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Abstract 



The study of anyons offers one of the most exciting challenges in contemporary physics. 
Anyons are exotic quasiparticles with non-trivial exchange statistics, which makes them 
difficult to simulate. However, they are of great interest as some species offer the prospect 
of a highly fault-tolerant form of universal quantum computation, and it has been sug- 
gested that the simplest such species may appear in the fractional quantum Hall state 
with filling fraction u — 12/5. Despite the current strong interest in the development of 
practical quantum computing, our ability to study the collective behaviour of systems of 
anyons remains limited. 

Meanwhile, tensor network algorithms are a relatively recent development in the field 
of condensed matter physics. They consist of Ansatze for the low-energy states of a 
lattice system, whose number of free parameters scale at most polynomially in the system 
size, and algorithms for their optimisation, manipulation, and analysis. However, many 
condensed matter systems possess a high degree of symmetry, which may be exploited 
to yield an even more efficient description of the low-energy subspace, and when I began 
work on my Thesis these algorithms (with the exception of DMRG) did not in general 
take advantage of these symmetries. 

In this Thesis I develop a formalism, based on the frameworks of spin networks and 
category theory, whereby a tensor network acting on any lattice model exhibiting a math- 
ematical structure corresponding to a Unitary Braided Tensor Category (UBTC) may be 
represented in a particularly compact and efficient manner corresponding to the exploita- 
tion of this structure. This permits the exploitation of global Abelian and non-Abelian 
internal group symmetries, both to facilitate the study of particular symmetry sectors of 
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the model, and for computational gain. Furthermore, the formalism also naturally ad- 
mits the study of models possessing non-trivial exchange statistics (e.g. fermions, Abelian 

anyons) and models possessing a UBTC structure which is not associated with a group 
(some Abelian and non- Abelian anyons), all for a computational cost polynomial in the 
system size. 

In addition 1 also describe the development of a tensor network algorithm to exploit the 
spatial symmetry of scale invariance present in quantum critical lattice systems. The 
resulting Ansatz provides a remarkably efficient description of the low-energy subspace of 
an infinite quantum critical lattice model, naturally yielding the polynomial correlators 
typical of such a system, and providing easy access to the majority of the conformal data 
which describe its behaviour in the continuum limit. Combining this Ansatz with the 
UBTC formalism for tensor networks provides a demonstration of the flexibility of these 
techniques, computing the conformal data associated with the continuum limit of two 
non- Abelian anyonic quantum critical lattice models. 

In summary, this Thesis provides a new Ansatz for the study of quantum critical lattice 
models, and a formalism permitting the exploitation of Abelian and non- Abelian symme- 
tries of lattice models, allowing the analysis of many fermionic and anyonic systems in 
polynomial time. 

No Ph.D. students were harmed during the making of this Thesis. 

Keywords: tensor networks, scale invariance, conformal field theory, symmetry groups, 
anyons 

Australian and New Zealand Standard Research Classifications (ANZSRC): 

020401 Condensed Matter Characterisation Technique Development (50%) 
020603 Quantum Information, Computation and Communication (50%). 
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Chapter 1 

Introduction 



1.1 How to Read This Thesis 

This Thesis presents a number of recent developments in the formahsm of tensor net- 
work Ansatze and algorithms, in which the author played a leading role. Tensor network 
Ansatze (or simply "tensor networks") are mathematical tools which may be used to 
efficiently represent a portion of the Hilbert space of a quantum mechanical system, fre- 
quently either the ground state or the low-energy subspace with respect to a specified 
Hamiltonian, and tensor network algorithms are algorithms for the efficient construction 
and manipulation of these tensor networks. They are frequently used to calculate proper- 
ties such as the low-energy spectra of lattice Hamiltonians, and the evolution of states as 
a function of time (see Sees. 1.3 and 1.4 for introductory citations). However, for a given 
Hamiltonian the ability of a tensor network algorithm to accurately represent the low 
energy subspace depends upon a number of factors. Existing tensor network algorithms 
tend to favour Hamiltonians which may be written as local operators over a small number 
of adjacent sites on a ID or 2D lattice, with other factors which determine whether a 
particular system may be efficiently analysed including the statistics of the system (e.g. 
bosonic, fermionic), and whether the structure of the tensor network reflects the structure 
of entanglement in the states being studied (Evenbly and Vidal, in preparation). 

The developments presented in the following chapters greatly extend the range of systems 
to which tensor network techniques may be applied, including formalisms for the study of 
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infinite systems (Chapter 2), fermions (Chapter 3), and anyons (Chapter 4). They also 
expand the capabihties of existing tensor network Ansatze when analysing symmetric 
Hamiltonians (Chapters 3 and 5), and yield substantial improvements in computational 
performance (Chapters 2 — 5). 

Although continuous tensor network algorithms do exist (Vcrstracte and Cirac, 2010), in 
this Thesis I will consider only tensor networks for lattice models. With the exception of 
Chapter 2, however, most of the results presented in this Thesis are completely general 
and may in principle be applied to any tensor network Ansatz or algorithm. However, 
when providing examples and demonstrations I will favour tensor networks of the Multi- 
scale Entanglement Renormalisation Ansatz (MERA) type. I therefore begin with a brief 
review of tensor network Ansatze in general and the ID MERA in particular. 

As described in the Statements of Contributions, much of the work presented in this Thesis 
has previously been published in international peer-reviewed journals. Rather than re- 
invent the wheel, this material has largely been reproduced verbatim in the Thesis. As 
a result the Thesis has a modular structure, where individual Chapters of the Thesis 
are essentially self-contained. Inevitably this approach comes at the expense of some 
repetition of background material, with individual Chapters, and sometimes Sections 
within those Chapters, often having their own introduction drawing the reader's attention 
to the relevant parts of this material. These are supplemented by material in the present 
Chapter, which provides context for the Thesis as a whole. Supplementary text at the 
beginnings or ends of the Chapters serves to bring out connections between the different 
research topics, and to place the individual research areas into a larger context. 

A Note on the Use of Personal Pronouns in this Thesis 

In this Thesis, use has been made of both the singular personal pronoun ("I"), and the 
plural ("we"), with the latter being used in different contexts to indicate either the author 
and the reader, or the author and his collaborators. Choice of personal pronoun should 
not therefore be treated as an indicator of whether work was performed independently or 
in collaboration. This information may be found in the Statements of Contributions in 
the preface to this text. 
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1.2 A General Introduction to Tensor Networks 

The idea that the state of a lattice model may be represented by a network of tensors may 
be motivated as follows: First, consider an n-site lattice Cq. If each site of this lattice 
is described by a rf-dimensional Hilbert space Hsite, then the Hilbert space of the n-site 
lattice is Hiattice = ('Hsitc)®"'- We may write a general state on lattice Cq as 

d 

1^) = 5Z (1.1) 

n,...,jn=l 

where \1) . . .\d) constitute an orthonormal basis of T^site, and |zi, . . . ,in) = Ki)® ^2)®- • 
\in)- For a fixed basis of "Hsite, a state of the lattice may be entirely specified by giving the 
tensor Q^..i„. It is convenient to introduce a graphical notation whereby a tensor with x 
indices may be represented by a blob with x legs, so that Q^..^^ is represented graphically 
by 




il 12 13 ■■■ in 



Optionally we may choose to give each leg a vertical orientation and identify upgoing 
legs with upper indices and downgoing legs with lower indices, but this practice is seldom 
followed for tensor networks which do not exploit internal symmetries as the metric, and 
hence conversion between upgoing and downgoing legs, is trivial. This notation may 
be understood as a simplified form of the Penrose graphical calculus (Penrose, 1971a). 
The distinction between upgoing and downgoing legs will become more important in 
the scheme presented in Chapters 4-5 for systems with non-Abelian symmetries and for 
anyons, and a different but related graphical notation will be introduced in Chapter 4 for 
the description of tensor networks with non-Abelian symmetries. 

The product of multiple tensors may similarly be represented by multiple blobs, with 
open legs corresponding to free indices and shared legs corresponding to summed indices. 
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For example, matrix multiplication may be represented as 

A^B^,^ Q (b^ _ (1.3) 

Unless otherwise specified, in this Thesis I will employ the Einstein summation convention, 
so /3 is summed over in the left-hand side of the above expression. If we now evaluate this 
sum over /3, we can write AgB^ = C", 




The idea behind tensor network Ansatze is that for any tensor Cj^ . j^ describing a state on 
a lattice, we may write down an equivalent collection of tensors linked by summed indices, 
such that when all these sums are evaluated, we recover the original tensor Cij...j„. Indeed, 
Eq. (1.4) may be thought of as a simple example where A and B form a tensor network 
which evaluates to give C. In general a tensor network diagram will contain multiple 
separate tensors, and the process of evaluating this diagram to obtain a single tensor 
(or, where there are no free indices, a number) is known as tensor network contraction. 
Similarly, taking two tensors within such a network and replacing them by a single, 
equivalent tensor (such that on contracting the entire network, the same resulting tensor 
is obtained) is termed contracting these two tensors. For example, A and B in Eq. (1.4) 
might in fact constitute part of a larger tensor network, such as the simple one shown 
below, and Eq. (1.4) then describes the contraction of A with B such that 




(1.5) 

In general, any tensor network may be completely contracted by a series of pairwise 
contractions of its component tensors. 

But why use tensor networks at all? The answer is simple: Efficiency. If a tensor network 
is to be capable of representing any state in the Hilbert space "H lattice, then it must 
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contain at least as many free parameters as Cjj...j„, i.e. ci", and in general will be even less 
convenient for computation than the form of Eq. (1.1). However, there exist particular 
choices of tensor network having less than c?" parameters which are nevertheless capable 
of providing an accurate description of an interesting subregion of this Hilbert space, 
for example the low energy subspace with respect to a particular Hamiltonian. Using 
these tensor networks, we may therefore numerically study the properties of a physical 
system governed by this Hamiltonian, typically at a fraction of the computational cost we 
would have incurred if we had chosen to retain the full description afforded by Cj^ . j^ in 
Eq. (1.1). This Thesis assumes a basic familiarity with the use of tensor network states 
and their associated algorithms, although for the reader desiring further material, a brief 
recapitulation of the MERA is provided in Sec. 1.3, and a number of references for further 
reading are listed in Sees. 1.3 and 1.4. 



1.3 Multi-scale Entanglement Renormalisation Ansatze 

Introduced in Vidal (2007a) and Vidal (2008), the family of tensor network Ansatze known 
as MERA are motivated by the idea of implementing a real-space renormalisation group 
transformation on the lattice. They represent a state using a layered structure, where 
each layer i may be considered as map between an initial ni_i-site lattice and a 
coarse-grained nj-site lattice £j, for n^-i > rii. In addition to the tensors which perform 
this coarse-graining, each layer of a MERA also incorporates a number of unitary tensors 
which act on the lattice to remove entanglement before coarse-graining takes place. 

As an example, consider the 3:1 ID MERA (Fig. 1.1). The open bonds at the bottom 
of Fig. l.l(i) correspond to indices ii, . . . ,ii8, the physical sites of the lattice. The open 
index at the top of the diagram enables the MERA to represent multiple sites within the 
Hilbert space, where the value of this index enumerates the represented states. This index 
ranges from 1 to Xtop- For a fixed value of this index the network represents a single 
state IV''') and may be contracted to the tensor Ci^...iis which specifies the coefficients of 
this state as per Eq. (1.1). A MERA having Xtop = 1 consequently represents only a single 
state Finally, we limit the dimension of each of the summed indices to at most Xi 
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(i) 




disentanglers 

Co 



(ii) 



(iii) 




Figure 1.1 - (i) Tensor network for a 3:1 ID MERA representing a subspace of dimension 
Xtop on a lattice of 18 sites. The top index of the tensor network ranges from 1 to xtop- 
(ii) Constraints on the disentanglers and isometry tensors of the 3:1 ID MERA, denoted u 
and w respectively, (iii) A reduced density matrix p may be constructed from the top tensor 
T and its conjugate as shown. 
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where x is a tunable parameter determining the number of free parameters in the Ansatz. 

The diagrammatic counterpart of Hermitian conjugation is implemented by vertically 
reflecting a tensor and complex conjugating all of its entries. Thus if we denote a dis- 
entangler by u"^ and its Hermitian conjugate by v}^^ = (m"^)*, then their diagrammatic 
representations are 

a 13 7 (5 



(1.6) 

7 5 a [3 

respectively. Similarly, the MERA for a bra, is constructed by vertically reflecting 
the ket, of Fig. l.l(i), and complex conjugating the coefficients of all the tensors in 
the network. 

For discussion of MERA algorithms to approximate the ground state of a system, see 
Dawson ct al. (2008), Rizzi ct al. (2008), and Evenbly and Vidal (2009a). A pedagogical 
introduction to the MERA, predominantly in one dimension, may be found in Vidal 
(2010), presenting interpretations of the MERA formalism both in terms of the real- 
space renormalisation group transformation (coarse-graining) described in this section, 
and also as a quantum circuit. Applications of the MERA to 2D systems may be found in 
e.g. Evenbly and Vidal (2009a), Evenbly and Vidal (2009b), Evenbly and Vidal (2010a), 
Evenbly and Vidal (2010b), Evenbly and Vidal (2010c), Cincio et al. (2008), Aguado and 
Vidal (2008), and Konig ct al. (2009). Note that Vidal (2010) also includes material on 
the scale-invariant MERA, which is the subject of Chapter 2 of this Thesis. 

1.4 Other Tensor Networks 

There also exist a number of other tensor network Ansatze and algorithms. This section 
lists a selection of introductory references and example papers for a few of the more 
popular. 



u 



and 
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1.4.1 Matrix Product States 

One of the most common tensor network Ansatze in use today is the Matrix Product 
State (MPS), which is the Ansatz underlying the Density Matrix Renormahsation Group 
(DMRG) technique developed by White (1992) for computation of ground states. Time 
evolution may be simulated using the Time Evolving Block Decimation (TEBD) algo- 
rithm of Vidal (2004). For further reading on DMRG and MPS, see Noack and White 
(1993); Perez-Garcia et al. (2007); SchoUwock (2005a, 2011); White and Noack (1992); 
White (1992, 1993); White and Feiguin (2004), and for TEBD of infinite chains, see Vidal 
(2007a). 

1.4.2 Tree Tensor Networks 

The Tree Tensor Network (TTN) may be used to represent states on lattices of arbitrary 
dimension. It has a heirarchical structure, but is not well suited to the representation of 
large critical systems due to the need for large bond dimensions (indices with large ranges) 
towards the top of the tree. Structurally, a TTN may be thought of as a MERA without 
disentanglers (Fig. 1.2). For further reading, see Shi et al. (2006), and also Fannes et al. 
(1992b); Friedman (1997); Lepetit et al. (2000); Martm-Delgado et al. (2002); Nagaj et al. 
(2008); Niggemann et al. (1997); Otsuka (1996); Taghacozzo et al. (2009). 

1.4.3 Projected Entangled Pair States 

The Projected Entangled Pair State (PEPS) Ansatz is a generalisation of the Matrix 
Product State Ansatz to two dimensions (Verstraete and Cirac, 2004). See also Jordan 
et al. (2008); Maeshima et al. (2001); Murg et al. (2007, 2009); Nishio et al. (2004); Sierra 
and Martm-Delgado (1998). 
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(ii) 




isometries 



Figure 1.2 - (i) Tensor network for a 3:1 ID TTN representing a subspace of dimension xtop 
on a lattice of 18 sites. The top index ( of the tensor network ranges from 1 to Xtop- Compare 
with the ID MERA of Fig. 1.1. (ii) Constraint on the isometries of the TTN (denoted w). 



10 Introduction 



Chapter 2 



Scale-Invariant MERA 



2.1 Introduction 

In this Chapter of the Thesis, we see how MERA-type tensor network Ansatze may 
be used to study the properties of infinite scale-invariant systems. Right from the start 
(Vidal, 2007a), the MERA has been constructed to implement a real-space renormalisation 
group (RSRG) transformation. It is by exploiting this property of the MERA that an 
Ansatz for scale-invariant systems may be constructed, as described by Giovannetti et al. 
(2008) and Pfeifer et al. (2009). Further publications studying and applying the scale- 
invariant MERA include those by Evenbly and Vidal (2009b), Montangero et al. (2009), 
and Giovannetti et al. (2009). Evenbly et al. (2010) address the application of the scale- 
invariant MERA to half-infinite and bounded ID chains, and most recently, Pfeifer et al. 
(2010) applies the scale- invariant MERA to a quantum critical system of anyons (see also 
Chapter 4 of this Thesis). 

The material presented in this Chapter is based upon research first published as Pfeifer, 
Evenbly, and Vidal (2009). The numerical results of Sec. 2.3, including Table 2.1 and 
Fig. 2.11, are reproduced or adapted from this reference and are (c) (2009) by the American 
Physical Society. 
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2.1.1 Real-Space Renormalisation Group Transformations 

Real-space renormalisation group transformations have a long history in condensed matter 
physics, dating back to Kadanoff 's spin-blocking technique (Kadanoff, 1966) and Wilson's 
solution of the Kondo problem (Wilson, 1975). However, such techniques really came of 
age with the development of the DMRG algorithm by White (1992, 1993). The defining 
feature of such techniques is that there exists some procedure whereby a theory on an 
initial lattice Co may be subject to some numerical coarse-graining process to yield an 
effective description on a new lattice £i, where each lattice site i on Ci corresponds to 
some region on £o; say sites ji, . . . , jn, and 

n 

dim(z) < J]dim(j,). (2.1) 

a=l 

This concept was originally proposed by Kadanoff (1966) in a classical context, with the 
idea of replacing a group of spins with a single effective spin chosen to be representative 
of the group. The first successful quantum mechanical application of an RSRG approach 
was the treatment of the Kondo problem by Wilson (1975), in which a coarse-graining 
procedure was chosen so that the retained portion of the Hilbert space corresponded 
locally to the low-energy eigenstates of individual terms of the Hamiltonian, e.g. in 
a Hamiltonian of the form 

H = J2k^+l (2-2) 

i 

where i and i+1 are sites of a ID lattice. However, the approach did not appear to 
generalise well to other problems. 

The development of DMRG by White in 1992 provided the next crucial insight — that 
the ground state wavefunction which minimises (H) does not necessarily also minimise 
the expectation value of each local term (/ij j+i), and the retained portion of the Hilbert 
space must therefore be chosen in a way which takes into account the total Hamiltonian. 
DMRG is one algorithm which satisfies this requirement. 

More generally, any tensor network which admits a description as a procedure mapping 
between a series of increasingly coarse-grained infinite lattices may be understood as 
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defining an RSRG transformation, and one can define a cost function [e.g. Tr^pH)] 
whose minimum corresponds to the desired state or subspace, and attempt to numerically 
optimise the defined RSRG transformation so as to extremise this cost function. For 
example, applying this philosophy to construct a quantum mechanical version of the 
spin-blocking technique, one obtains the Tree Tensor Network formalism of Shi et al. 
(2006) (see Sec. 1.4.2, above). 

When an RSRG transformation is applied to a Hamiltonian Hq on a lattice Co, it yields 
a Hamiltonian Hi on the coarse-grained lattice Ci. Repeated application of the RSRG 
transformation therefore causes the set of Hamiltonians Hi to describe a trajectory in the 
space of Hamiltonians, termed a Renormalisation Group (RG) flow. If the Hamiltonians 
Hq and Hi satisfy Hi = Hq, then the Hamiltonian remains unchanged under repeated 
application of the RSRG transformation, and we term Hq a fixed point of the RG flow. 

2.1.2 Lattice Models Exhibiting Scale Invariance 

In order to exhibit scale invariance, a lattice model must be free of any characteristic 
length scales. Consequently, when a Hamiltonian is a fixed point of the RG flow defined 
by an RSRG transformation TZ, then the correlation length ^ for all operators in a system 
must be either zero or infinite. 

Recall now that our tensor network Ansatz must not only describe the RSRG transforma- 
tion 71, but also provide an accurate description of the low-energy subspace of the Hilbert 
space of the system. Let us consider the different possible types of scale-invariant lattice 
model we may encounter: 

1. ^ = 0, ground state is a product state: Unentangled. Product states may trivially 
be described by any tensor network. 

2. ^ = 0, ground state is a topologically ordered state: Entangled. Experience indicates 
that such states may be efficiently described by a MERA with finite bond dimension 
(Aguado and Vidal, 2008; K5nig et al., 2009). 

3. ^ = oo, quantum critical system: Highly entangled. In the ground state, entangle- 
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ment of a contiguous region A of length L with the rest of the lattice (as measured by 
the von Neumann entropy S) scales as 5* ~ log(L) in ID, and S* ~ L or 5* ~ L log(L) 
in 2D, depending on the model under consideration. 

Obviously, we may simulate a product state using any tensor network we like, including 
the MERA. Prior experience shows that for topologically ordered systems, the MERA 
may once again be a good choice. Finally, what about quantum critical systems? We 
require a tensor network capable of encoding a bipartite entanglement entropy which in 
ID should scale as log(L) when evaluated for a contiguous region A^^^ having length L, 
and in 2D should scale at least as L, and preferably as Llog(L), for a region A'^'^^ having 
dimension L x L. To understand why the MERA is once again the natural choice, it is 
necessary to briefly examine how entanglement entropy scales for any tensor network. 

Consider a state \ip) represented by a single tensor Q^ . j^ as in Eq. (1.1). Suppose we 
wish to investigate the entanglement between two regions of the lattice, sites ii . . .ia and 
ia+i . . - in- If we perform a Schmidt decomposition of state then we write 



where \j is a list of strictly positive coefficients. In terms of linear algebra, this is 
equivalent to performing a singular value decomposition on a matrix Ck^k2-, where in- 
dices ki and k2 enumerate states on indices ii, . . . ,ia and ia+i, ■ ■ ■ ,in respectively. The 
dimension of index j may range from 1 (for a product state) to the lesser of dim(fci) = 
dim(zi) X ... X dim(za) and dim(fc2) = dim{ia+i) . . . x dim(z„) (for a highly entangled 
state). Assuming that the coefficients Xj are sorted in decreasing order of magnitude, 
then the more entangled the state, the higher the value of j before Xj ^ Ai. 

Now, consider as an example a ID state represented by an MPS, where the range of the 
indices in the network has been limited to x- If we perform a bipartition of such a state, 
as shown in the Fig. 2.1(i), then the range of j in Eq. (2.3) will be limited to at most X; 
and this provides an upper bound on the amount of entanglement which such a state may 
represent. 




(2.3) 



In general, a crude quantification of the maximum entanglement a tensor network may 
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encode between a region A and the rest of the lattice is therefore given by taking the 
product of the dimensions of the bonds which one must cut to separate the tensor network 
into two regions, one contacting the physical lattice only within region A, and the other 
contacting the lattice only outside of A [e.g. Fig. 2.1(ii)]. (Taking the logarithm of one 
over this value yields an upper bound on the von Neumann or entanglement entropy, S.) 
In general, multiple such cuts exist, and the maximum amount of entanglement which 
may be encoded is determined by the cut giving the smallest value. 

By studying the scaling of S with the size of region A, we observe that for an MPS, the 
maximum amount of entanglement which may be encoded is independent of the size of 
region A, S ^ const. For a TTN, the amount of entanglement which may be encoded 
exhibits a more complicated dependence on the size and position of region A. However, 
for an n-into-1 tree it is possible for any k G Z+ to choose a region A of linear size 
L = kn, having the same entropy as an appropriately chosen region of length L = n. The 
overall performance of the Ansatz is limited by this worst-case scenario, and consequently 
a TTN of constant x ^i^lso exhibits an entropy scaling S ~ const. In contrast, the ID 
MERA exhibits a scaling S ~ log(L), making it well-suited to the study of quantum 
critical systems. In 2D the situation is less ideal, with the 2D MERA exhibiting a scaling 
5* ~ L, and thus only being suited to the study of critical systems which do not display 
a logarithmic correction to the entanglement entropy. However, recent development of a 
branching MERA algorithm (Evenbly and Vidal, in preparation) with entropy scaling as 
S ^ L log(L) or better suggests that once again, some form of the MERA may prove to 
be a good choice of Ansatz for all critical systems in 2D. 

This is not to say that the MPS and TTN cannot be used to calculate properties of quan- 
tum critical systems. They can, and in many situations may yield excellent numerical 
approximations to ground state energies and short-range correlators. However, due to 
their limited capacity for encoding entanglement within the structure of the tensor net- 
work, attempts to construct a RSRG transformation will fail at sufficiently large length 
scales (and we shall see in Sec. 2.2.2 that many interesting properties of quantum critical 
systems may be computed in the large- length-scale, or infra-red, limit). Consider as an 
example the application of an n-into-1 TTN to a highly entangled lattice model on a 
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(i) 



(ii) 



(iii) 



-III— 

cut 



B 



ooo 




ooo 



cut 
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B2 



Bo 
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Figure 2.1 - (i) Bipartition of an MPS into two regions A and B. The dimension of the bond 
traversing the partitioning of the system limits the amount of entanglement between the two 
regions, (ii) Entanglement of a region A with the rest of the system (5i and B2) is limited 
by the product of the dimensions of the bonds which must be cut to separate A from the 
rest of the system. Multiple such cuts exist for a given region A, for example the alternative 
shown in (iii), each placing an upper bound on the amount of entanglement. Only the most 
stringent such bound is therefore of interest. 
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lattice Cq. A single layer of the TTN coarse-grains Cq into an effective lattice £i, and 
each site in Ci now corresponds directly to n sites in Cq. The entanglement between a 
single site of Ci and the rest of the lattice will thus be the same as between those n sites 
on £o and the rest of the lattice. If the entanglement entropy of the model exhibits any 
dependence on L, then repeated coarse-graining of the lattice will cause this entangle- 
ment to continually increase. If the dimension of the indices of the TTN is bounded by 
some value Xi then after some number of coarse-graining steps, the approximation made 
in imposing this limit on index dimension will lead to a failure of the TTN to accurately 
reproduce the properties of the ground state over large length scales. Alternatively, the 
index dimension would have to increase with each layer of the coarse-graining procedure, 
eventually becoming infinite. First, this is computationally unfeasible, and second, we 
anyway desire that the structure of our tensor network should reflect the scale-invariant 
nature of the ground state. 

In contrast, the MERA may be thought of as a TTN supplemented by additional tensors 
known as disentanglers (for illustration of this in ID, compare Figs. 1.1 and 1.2). When 
the MERA is interpreted as an RSRG transformation, then the disentanglers in each level 
act to remove short-range entanglement from the ground state. We anticipated that in 
ID, for a well-optimised MERA representation of the low-energy subspace of a quantum 
critical system, they would do so to a sufficient extent that the entanglement entropy of 
a region A' of length L on lattice would be the same as the entanglement entropy of 
a region A of length L on lattice Cq, for any coarse-grained lattice C^, even when the 
entanglement entropy on an individual lattice scales as ~ log(L). This proposition was 
based on the observation that the ID MERA is constructed to be capable of encoding an 
entanglement entropy which scales as S* ~ log(L), and its predictions have been borne 
out by subsequent experience. 



2.1.3 Interesting Properties of ID Quantum Critical Systems 

In studying quantum critical systems, we are particularly interested in their behaviours in 
the infra-red limit. When we take the continuum limit of a ID quantum critical system, 
we obtain a Conformal Field Theory (CFT) in 1+lD (Cardy, 1996; Di Francesco et al.. 
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1997) which describes the infra-red behaviour of the system in the vicinity of the associated 
phase transition, and a useful question to ask is whether we can extract from our Ansatz 
sufficient data to identify and fully characterise the CFT. 

In the operator formalism, a 1+lD CFT may be described in terms of an infinite number 
of operator- valued fields (l)a{x,t). It is conventional to first define the theory on the 
cylinder, with t G [—00, +00] and x G [0, L), before mapping to the complex plane via the 
reparameterisation 

with the fields of the theory now being denoted (pai^, z). Under the action of a conformal 
mapping z ^ z' = f{z), the correlators of these fields transform as 

z,)M^2, Z2)...) = l[ f{z^)^"f{z.t-{M^[, z[)M4, 4) • • •), (2.5) 

a 

where each field (pa is associated with a holomorphic and an antiholomorphic conformal 
dimension, ha and ha respectively. These in turn may be combined to give the scaling 
dimension of the field, Aq, = ha + ha, and the conformal spin, Sa = ha — ha- 

For any 1+lD CFT, these fields (pa, which we will term scaling fields, may be organised into 
conformal families, each consisting of an infinite number of fields. Within each conformal 
family, the field with the smallest scaling dimension is termed the primary field, with all 
others being termed descendants. We may associate with each operator field (pa a state 
I0Q,) generated by acting with (pa on the origin of the vacuum state (which corresponds 
to t = — oo), 

10.) =0a(O,O)|O). (2.6) 

If we define the operators 

Ln = ^j z^^'f{z)dz, Ln = ^j z^^^f{-z)d-z, (2.7) 

where T is the energy-momentum tensor, and the contour integration is performed over 
any contour which encircles the origin, then these operators I/„ and L„ form representa- 



2.1 Introduction 



19 



tions of the Virasoro algebra. They obey the commutation relations 

[Lm, Ln] = (m - n)Lm+„ + —{w? - m)5m+n,o, (2.8) 
[Lm, L„] = (m - n)L„i+n + — (m^ - m)5„+„,o, (2.9) 
[L^,L„] = 0, (2.10) 

where parameter c is a constant known as the central charge of the CFT, and repeated 
application of L„, n < and Lm, m < to the state |0q,) associated with the primary 
field of any conformal family will generate all other states associated with members of 
that family. There may be a finite or an infinite number of conformal families, but of 
greatest interest to us will be the CFTs known as minimal models, for which the number 
of conformal families is finite. 

The identity operator is always the primary field for one of the conformal families, and 
we may always choose our operator fields to satisfy 

{4>a{z,z))=S^,, (2.11) 

C 

where Cap = 5aj3 corresponds to a particular choice of normalisation, and 5ai = 1 if 0q, = I 
and otherwise. We must also specify the coefficients C^^^ of the three-point function. 



ha-\-hi^ — hp-\-h'y — ha hj-l-ha — hp 



— hc(.-\-h^ — /i^ _ h0-\-h'y — /iq: _ h^-{-h(x — 
Z^p Z^^ Z^a 

Zap = \Za — Zp\ Zap = \Za — Zp\. (2-14) 



Whereas Cap was a normalisation factor which we were free to choose as we liked, the 
values of Cap-^ form part of the description of the CFT under consideration. 

Expression (2.13) implies an algebra known as the Operator Product Expansion (OPE), 

^a{Za,Za)^p{zp,Zp) = ^C aP^{Za - Zp)'^'^-^^^^'' {Za - ZpY^'^'^^^^'' ^^^{Za,Za) + . . . , 

(2.15) 
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which may be inserted into higher-order correlation functions with higher-order terms 
vanishing in the hmit that \za — Zj^l is much smaller than any other separation in the 
correlator. 

To fully describe a 1-l-lD CFT in the operator formalism, it suffices to specify 

1. The primary fields, (pa- 

2. Their scaling dimensions A^^ and conformal spins Sq,. 

3. The central charge c of the Virasoro algebra. 

4. The coefficients Ca/sj of the operator algebra for the primary fields. 

We will see that it is possible to extract all of these data from the Scale- Invariant MERA, 
with the exception of the conformal spin. However, the data which can be obtained are 
nevertheless frequently sufficient to uniquely identify the CFT describing the infra-red 
limit of a particular quantum critical system. 

2.2 Scale- Invariant MERA Algorithm 

In Sec. 2.2.1, I describe the algorithm for constructing an infinite, scale- invariant MERA. 
This approach may be applied to either the ID or the 2D MERA, but in this Thesis I 
will primarily address the study of ID quantum critical systems, whose infra-red limits 
correspond to the interesting and highly-studied I+ID CFTs. Material on the extraction 
of conformal data in Sec. 2.2.2 is addressed primarily to these systems, and to computation 
of the parameters described in Sec. 2.1.3. 

An example application of the Scale-Invariant MERA algorithm to infinite 2D lattice 
models may be found in Evenbly and Vidal (2009b). 
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2.2.1 Construction of MERA for the Low Energy Subspace 
2.2.1.1 Overview 

For a finite system, a MERA normally consists of a finite number of layers of tensors, 
each layer consisting of a row of disentanglers and a row of isometries (Fig. 1.1). Each 
layer performs a coarse-graining procedure, mapping from a lattice to a coarser 
lattice £j. This process incorporates a truncation of the Hilbert space, such that after 
all layers of the MERA have been applied, the dimension of the Hilbert space on the 
maximally coarse-grained lattice is sufficiently small to exactly diagonalise. Numerical 
optimisation of the MERA (Dawson et al., 2008; Evenbly and Vidal, 2009a; Rizzi et al., 
2008) is performed to ensure that the Hilbert space of the final coarse-grained lattice 
exhibits maximal overlap with the interesting region of the Hilbert space of the original 
lattice, typically the low-energy subspace of a system. 

For an infinite system, this procedure obviously requires some modification. No matter 
how many times we apply a coarse-graining transformation to an infinite lattice Cq, the 
result is always an infinite lattice, and the system never becomes small enough to exactly 
diagonalise. However, the scale-invariant property of quantum critical systems comes to 
our rescue. To see how this works, let us assume that we have a Hamiltonian Ho^ip which 
is constructed on lattice Cq and lies exactly at the fixed point of an RG fiow. 

We know that if we were able to construct a MERA with an infinite number of layers which 
represented the low-energy subspace of this Hamiltonian, then because -f^o.fp is a fixed 
point of the RG fiow, application of a layer of the MERA to perform a coarse-graining 
from Cq to Ci would map i^o.fp into an identical operator i^i,fp on the coarse-grained 
lattice. An object which maps operators into operators is termed a superoperator, and 
we may therefore define the scaling superoperator S as the superoperator implemented 
by this layer of the MERA, which maps operators from lattice Cq to Ci for our scale- 
invariant system. Of course, because Ho,ip and -ffi,fp are identical, the layer of MERA 
constructed on Ci will be identical to that constructed on Cq, and we may equally well 
define S with reference to any layer of this infinite MERA. Because the Hamiltonians are 
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similarly identical, we will drop the lattice index, and simply write H{p for the fixed-point 
Hamiltonian on any lattice £j. 

Because the Hamiltonian is identical on all lattices Ci, the reduced density matrix pfp 
which minimises the energy Tr(/3fpiffp) is similarly also identical on every layer of coarse- 
graining. However, in a MERA we may always calculate the reduced density matrix on 
a lattice from the reduced density matrix on lattice Ci (Evcnbly and Vidal, 2009a). 
Let us denote by S* the superoperator which is the dual of S, and maps operators on 
into operators on Because pfp is identical on every layer, it must be an eigenoperator 
of S*, and because Tr(pfp) = 1 on every layer, it must have eigenvalue 1. Provided there 
exists only one eigenoperator of S* which has eigenvalue 1, this then suffices to uniquely 
define the fixed-point reduced density matrix pfp. When S* has only one eigenoperator 
with eigenvalue 1, knowledge of the scaling superoperator S and its dual are sufficient to 
compute the reduced density matrix, and these superoperators in turn may be constructed 
from any layer of this infinite, scale-invariant MERA. 

Finally, because all layers of this MERA are identical, we need only describe the disen- 
tanglers and isometrics of one layer in order to describe the state of the entire system. 
Assuming also translation invariance, we need only one disentangler and one isometry in 
order to describe the entirety of this infinite MERA, or compute the reduced density ma- 
trix on any lattice Ci. What will be presented in this Section is therefore an algorithm for 
determining exactly these tensors: The disentangler and isometry of the scale-invariant 
MERA. 

2.2.1.2 A Less Idealised Situation 

In the above discussion, it was assumed that the Hamiltonian of the system was precisely 
the fixed point Hamiltonian of the RSRG transform; that is. Hi = Hq = Hfp. For this to 
be true, Hfp must correspond to a scaling field, or sum of scaling fields, of the associated 
CFT, all with identical scaling dimension A^. 

In practice, the Hamiltonian of the quantum critical system may not be exactly the fixed 
point Hamiltonian, but may also include additional scaling fields, provided these fields 
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have scaling dimension > A^. On repeated coarse-graining these fields are suppressed 
relative to the Hamiltonian. While these fields will in theory never vanish completely, 
and on repeated coarse-graining Hi will only approach the fixed point of the RG flow 
asymptotically, 

lim i^, ^fp, (2.16) 

we will assume that they decay sufficiently rapidly that their existence may be neglected 
after some finite number r of applications of the coarse-graining process. We therefore 
construct our Ansatz to consist of r layers of ordinary MERA, acting on lattices Cq to 
Cr-ii after which the difference between Ht on lattice Cr and H{p is negligible, and all 
subsequent layers of the MERA will be essentially identical. We therefore surmount the 
layers 1 . . . r of the MERA by one further layer t + 1, which acts on lattice and is 
assumed to be repeated an inffnite number of times (as layers r + 2 and above, acting 
on lattices £t+i to £oo)- It is this layer r + 1 of the MERA which is then used in the 
construction of the scaling superoperator. 

We will call an operator a scaling operator if it is an eigenoperator of S, and in Sees. 2.2.2- 
2.3 we will endeavour to identify these operators with the scaling ffelds 0q, of the CFT 
associated with our quantum critical lattice model. As a note of terminology, scaling fields 
which decay more rapidly than the fixed point Hamiltonian under the action of the RSRG 
transformation are termed irrelevant, as are the associated fields of CFT. Those which 
decay at the same rate are termed marginal, and those which decay less rapidly are termed 
relevant. We will adopt the same terminology for scaling operators. The Hamiltonian of 
a quantum critical system will only ever contain marginal and possibly irrelevant terms. 

2.2.1.3 The Algorithm Itself 

I now present explicitly a practical algorithm for optimisation of the Scale-Invariant 
MERA for a local quantum critical Hamiltonian on a lattice, which may contain irrelevant 
terms. This algorithm will be described in general language applicable to both ID and 
2D systems, though accompanying illustrations will refer specifically to the 3:1 MERA in 
ID. 

It is assumed that the Hamiltonian under consideration is nearest-neighbour, next-to- 
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nearest neighbour, etc., as appropriate to the MERA being employed (for example, the 3:1 
MERA on the ID lattice is constructed for the study of nearest-neighbour Hamiltonians) . 
This may always be achieved by means of some initial coarse-graining onto an effective 
lattice of higher site dimension if required. As an example, we will subsequently consider 
the critical Ising model on a ID lattice, which is a nearest-neighbour Hamiltonian and 
thus is directly suitable for analysis using the 3:1 MERA. 

The MERA is initially constructed to consist of some small number of free layers 1 . . . r, 
where r is typically around 1 or 2, over which the local dimension of the lattice increases 
from d on io x ^ d on £r-i- These are surmounted by the scale-invariant layer, 
which maps from lattice Cr-i to lattice C-r-, both of dimension x- Initial choices of tensor 
for the disentanglers and isometrics are comparatively unimportant, and may be chosen 
randomly within the constraints of Fig. l.l(ii), or assigned systematically to some known 
initial configuration. Optimisation then proceeds as follows: 

1. Construct the fixed-point reduced density matrix, Pfp = Pr+i on lattice Cr+i, by 
diagonalising the dual of the scaling superoperator S* (see Fig. 2.2 for an example). 

2. "Descend" p^+i to obtain . . . po in the usual manner (Evenbly and Vidal, 2009a). 

3. Proceeding row by row from layer 1 to layer r of the MERA, for each layer i: 

(a) Update the disentanglers in the usual manner (Evenbly and Vidal, 2009a). 

(b) Update the isometrics in the usual manner. 

(c) "Lift" the Hamiltonian from to in the usual manner. 

4. The Hamiltonian Hr+i now closely resembles the fixed point Hamiltonian Hf^. Op- 
tionally, we may now "lift" Hr+i a few more times, using the tensors of the scale- 
invariant layer r -|- 1 of the MERA, to obtain a tensor H*^i which is even more close 
to Hf^. Using H*j^i in lieu of Hr+i may yield more accurate computation of critical 
exponents, but at the cost of slower convergence of the Scale-Invariant MERA. 

5. Using Ht+i (or H*^i, if preferred) and the reduced density matrix p^+2 from one 
layer further up the infinite MERA (which is taken to be the same as Pr+i), update 
the disentanglers of layer r -|- 1 of the MERA (the scale-invariant layer). 
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6. Compute numerical properties (e.g. ground state energy, scaling dimensions; see 
Sec. 2.2.2). 



7. Repeat all steps until the chosen cost function is satisfactorily converged. 



Following initial convergence of the MERA, the quality of the numerical results (e.g. 
ground state energy, scaling dimensions, etc.) may be increased by adding more free 
layers below the scale-invariant layer. To do so, copy the tensors of the scale-invariant 
layer (denoted r + 1) to obtain a layer r + 2. Layer r + 2 is now the scale-invariant layer, 
and the above optimisation procedure is now repeated with layers 1 . . . r + 1 optimised 
in the usual manner for a standard MERA, and layer r + 2 being used to construct the 
fixed-point reduced density matrix and the scaling superoperator. This process may be 
repeated until insertion of additional layers no longer causes a significant change in the 
computed properties of the MERA. 

The above algorithm serves as an illustrative example as to how a scale-invariant MERA 
may be converged. In practice, a significant time saving may be made by modifying 
the computation of pr+i- Rather than computing the dominant eigenoperator of S* 
exactly on every iteration, we instead assume that after updating the scale-invariant 
MERA layer, pr+i from the previous iteration has a non-trivial overlap with the dominant 
eigenoperator of the new S* . We therefore take Pr+i from the previous iteration, and 
apply the dual of the scaling superoperator once (i.e. we descend this operator using 
the scale-invariant layer r + 1). We then take the resulting operator to be the new 
Pr+i ■ ■ - Poo- In the limit that the MERA converges (assuming, as always, that this limit 
exists — an assumption borne out well in practice), S* remains constant from iteration to 
iteration, and thus is repeatedly applied to Pr+i, which will thus gradually converge to the 
dominant eigenoperator of S* as required. In practice, this process leads to a co-ordinated 
convergence of S* and Pr+i, and requires less time than exactly computing pr+i on every 
iteration. 
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Figure 2.2 - Construction of the two-site fixed-point reduced density matrix pfp of a ID 3:1 
MERA by diagonalising the dual of the scaling superoperator S*. The reduced density matrix 
is found by solving the above graphical equation, where the disentanglers and isometries come 
from the scale-invariant layer r -|- 1 of the MERA. We assume that the solution is unique (i.e. 
that S* has precisely one eigenoperator with eigenvalue 1), an assumption which proves valid 
for the systems studied in this Thesis. A reduced density matrix which satisfies this equation 
remains unchanged upon being "descended" from layer to layer within the scale-invariant 
region of the MERA. 
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2.2.2 Extraction of Conformal Data 



As described in Sec. 2.1.3, we may associate the infra-red (large-scale) behaviour of a ID 
quantum critical theory with a I + ID CFT. To extract the conformal data describing this 
CFT, we must first identify the objects in the quantum critical theory which correspond 
to the scaling fields of the CFT. These are objects which remain invariant under the 
action of an RSRG transformation, and consequently may be identified with operators 
which are eigenoperators of the scaling superoperator. Note that in this Section, we are 
interested in calculating properties in the large-scale regime of the quantum system, and 
consequently all disentanglers, isometries, reduced density matrices, etc. are drawn from 
the scale-invariant layer of the MERA, which is assumed to be repeated an infinite number 
of times and therefore describes the behaviour of the quantum critical system on all larger 
length scales. 

For the 3:1 MERA, the causal cone has a width of two sites, and consequently we may 
construct a two-site scaling superoperator 5(2) [Fig. 2.3(i)] whose eigenoperators are two- 
site scaling operators [Fig. 2.3(ii)]. Note that the scaling superoperator is the average of 
three diagrams. This is because a two-site operator on the coarse-grained lattice receives 
contributions from operators on three distinct pairs of sites on the fine-grained lattice, 
and this must be taken into account in the construction of the scaling superoperator. The 
two-site reduced density matrix is similarly an eigenoperator of <S'^2)^ dual of 5(2), with 
eigenvalue 1, as shown in Fig. 2.2. However, we also note that on privileged sites of the 
3:1 MERA, it is also possible to consider one-site scaling operators which remain invariant 
under the action of the RSRG transformation. These operators are eigenoperators of the 
one-site scaling superoperator, iS(i), as shown in Fig. 2.4. 

Note that due to the contraints on the disentanglers and isometries [Fig. l.l(ii)], both 
scaling superoperators are necessarily unital, S(I) = I, so that the identity operator is 
always a scaling operator with eigenvalue Ai = 1, and contractive, \Xa\ < 1 V (Bratteli 
and Robinson, 1979). 

As the sites of the MERA are spacelike-separated, the one- and two-site scaling operators 
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Figure 2.3 - (i) The two-site scaling superoperator 5(2) of the 3:1 IVIERA. It is constructed 
from the isometries and disentanglers of the scale-invariant layer of the MERA. (ii) Two-site 
local scaling operators are eigenoperators of 5(2) ■ Note that the scaling superoperator is the 
average of three diagrams. This is because a two-site operator on the coarse-grained lattice 
receives contributions from three distinct pairs of sites on the fine-grained lattice, and this 
must be taken into account in the construction of the scaling superoperator. 
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(i) 



(ii) 



5(1) - 



Figure 2.4 - (i) The one-site scaling superoperator, constructed from the isometries of 
the scale-invariant region of the MERA. (ii) One-site scaling operators are eigenoperators of 
5(1). 
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where purely spatial co-ordinate, and 



(2.17) 
(2.18) 

(2.19) 



(2.20) 



For convenience, we shall now choose to work specifically with the one-site scaling oper- 
ators. We may normalise these scaling operators by imposing Eq. (2.18) with Cq,^ = Sap- 
When two scaling operators are placed on consecutive lattice sites, the correlator (2.18) 
reduces to 

(0„(O)0^(1)) = S^^, (2.21) 

where the expectation value is computed with respect to the two-site reduced density 
matrix on the scale-invariant portion of the MERA, which we will denote p^^y For the 
3:1 MERA, this is the same as the fixed-point reduced density matrix pfp calculated 
during optimisation. The local scaling operators (pa niust therefore be (ortho) normalised 
to satisfy 

TT[{k^4^p)P(^)] = Sap, (2.22) 
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Figure 2.5 - Diagrammatic expression for the normalisation of (0q,(O) 



as represented graphically in Fig. 2.5. 

We may now compute correlators for pairs of one-site local scaling operators, provided 
these operators are located at a separation such that each application of the MERA maps 
a one-site local scaling operator into a one-site local scaling operator (e.g. Fig. 2.6). We 
find that under these conditions, the correlators scale as 

{kM^p{xp)) = 5„;3(A„)2i°S3(r-.,). (2.23) 

Using the identity 

^logy ^ ylogx (2.24) 

(try taking the log of both sides), we see that the scaling dimensions of the primary fields 
may be computed according to 

A„ = -loggA^. (2.25) 
(More generally, for an n-into-1 MERA, the same argument yields = — log^A^.) 

We may also use correlators to calculate the OPE coefficients Ca/s-y- Although we do 
not have direct access to a translation-invariant scale-invariant reduced density matrix 
p^3^ in direct analogy to the two-site reduced density matrix p^2)) we may nevertheless 
easily compute three-point correlators (0q-(xo — 1) k{xo) (/)-y{xo + 1)) for certain privileged 
locations on the lattice. Assuming translation invariance of (0a(xo — 1) (pplxo) 4>^{xo + 1)) 
and symmetrising across the two diagrams given in Fig. 2.7 suffices to give us the values 

of Cal3-y- 

C^fs^ = 2^-+^"-^''(0„(O)0«(l)0«(2)). (2.26) 
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Figure 2.6 - Example of a longer- ranged correlator: Computation of (^a(xo) ^^(xo -I- 9)) 
for one-site scaling operators on privileged sites with causal cones of width 1. Note that all 
tensors, on both layers, correspond to those of the scale-invariant region of the MERA. 
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Figure 2.7 - Graphical expressions for the three-point correlator {(f>a{xo — l) 4>/3{xo) 0^(xo+l)) 
at two privileged locations on the ID lattice. All disentanglers, isometries, and reduced density 
matrices are those from the scale-invariant region of the MERA. 



If sufficient computational resources are available, a more rigorous evaluation of Ca/s-y may 
be achieved by computing p^g^ according to Fig. 2.8, and then determining (0q,(O)0q,(1)0q,(2)) 
using p^3^ according to Fig. 2.9. 

We may also compute the central charge, which is obtained from the von Neumann 
entropies associated with the one- and two-site reduced density matrices of the scale- 
invariant layer of the MERA: 

5(p) = -Tr(plog2/3) (2.27) 
c = 3[S{p^^)-S{pf,^)]. (2.28) 

(The one-site scale-invariant reduced density matrix may be obtained by symmetrising 
over the two ways of tracing out one site of p^2) — see Fig. 2.10.) 

How well may these calculations be expected to work in practice? It is important to 
recognise that for any CFT, there are always an infinite number of local scaling opera- 
tors. However, for any MERA with finite bond dimension x, there are only ever a finite 
number of eigenoperators of the scaling superoperator. We anticipated, and this is borne 



2.2 Scale-Invariant MERA Algorithm 



33 




Figure 2.8 - Expression to be solved for p^^y All disentanglers, isometries, and reduced 
density matrices are those from the scale-invariant region of the MERA. 
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Figure 2.9 - Correlator {(j)a{0)4'a{'^)4>a{'^)) evaluated using 
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Figure 2.10 - The one-site scale-invariant reduced density matrix p^^^ may be obtained by 
symmetrising over the two different ways of tracing out one site of the two-site scale-invariant 
reduced density matrix p^^^y 



out in practice, that by constructing a MERA which represents well the low-energy sub- 
space of the quantum critical theory, we would obtain to a high level of accuracy the 
conformal parameters associated with the scaling fields of lowest scaling dimension, but 
that this accuracy would decrease on going to larger scaling dimensions, with an inevitable 
truncation at some finite scaling dimension A^ax- The space of states = 0a(O)|O) 
associated with the local scaling operators of the MERA is therefore a finite-dimensional 
vector space, on which exists at best only an approximate, truncated representation of 
the Virasoro algebra describing the associated CFT (see Eqs. 2.8-2.10 of Sec. 2.1.3), be- 
coming exact in the limit x ^ oo. Nevertheless, we find that even for relatively modest 
X, it is frequently possible to construct a MERA which yields reasonable accuracies for 
the conformal data. 
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Field A^^T ^MERA j^elative error (%) 

Ising model 
10 

a 1/8 = 0.125 0.124997 0.002 

e 1 1.0001 0.01 



Potts model 



I 












0"! 


2/15 = 


0.13 


0.1339 


0.4 




2/15 = 


0.13 


0.1339 


0.4 


e 


4/5 = 


0.8 


0.8204 


2.5 




4/3 = 


1.3 


1.3346 


0.1 


Z2 


4/3 = 


1.3 


1.3351 


0.1 



Table 2.1 - Scaling dimensions, exact (A*-'^^) and computed (A'^^^^), for the primary fields 
of the ID Ising and three-state Potts models. Numerical results were obtained using a 3:1 
Scale-Invariant MERA with x = 22. Table adapted from Pfeifer, Evenbly, and Vidal, Physical 
Review A, 79, 040301, 2009, © (2009) by the American Physical Society. 



2.3 Results 

This Section presents results demonstrating the capabilities of the Scale- Invariant MERA. 
Two systems were studied: The Ising model, and the three-state Potts model, which are 
known to be associated with CFT minimal models A^(4, 3) and A^(6,5) respectively. 
Using a X = 22 MERA, scaling dimensions of the primary fields were obtained to within 
a relative error of 0.01% for the Ising model, and 2.5% for the three-state Potts model 
respectively (Table 2.1), with appropriate multiplicities for all primary fields and also the 
lower-scaling-dimension secondary fields (Fig. 2.11). 

The computed central charges of cising = 0.5007 and cpotts = 0.806 closely reflected the 
exact values of c£|^^ = 0.5 and Cp^^^g = 0.8 respectively. OPEs were also computed for the 
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Figure 2.11 - Scaling dimensions obtained from the spectra of the scaling superoperators 
S. Circles indicate primary fields. Left: For the Ising model we can identify the scaling 
dimensions of the three primary fields, the so-called identity (I), spin {a), and energy (e) 
fields, together with several of their descendants. Right: The spectrum of S for the three- 
state Potts model shows the primary fields I (An = 0), ai and cj2 {^ai = ^0-2 = 2/15), e 
(Ae = 4/5), and Zi and Z2 (A^^ = A^j = 4/3), along with the first descendants of cji and 
CJ2. Figure reproduced from Pfeifer, Evenbly, and Vidal, Physical Review A, 79, 040301, 2009, 
© (2009) by the American Physical Society. 




Ising model and compared with the exact figures of 

^CFT _ X riCFT _ 1 ^CFT _ ^CFT _ ^CFT _ n /'o on\ 

^al3I — Oap (-^acje " 2 ^'y'^^ ~ ^^^^ ~ ^'^^'^ ~ [^.Zy) 

and permutations of the indices thereon, with errors in all values being bounded by 3 x 10~^ 
(Pfeifer et al., 2009). 

These data are easily sufficient to identify the CFTs associated with these lattice models 
to a high degree of accuracy, and confirm the hypothesis that the Scale-Invariant MERA 
is an effective Ansatz for the description of quantum critical systems in one dimension. An 
example of the application of the scale-invariant MERA algorithm to infinite 2D lattice 
models may be found in Evenbly and Vidal (2009b). 



2.4 Adding a Boundary 



37 



2.4 Adding a Boundary 



The author of this Thesis was also briefly involved in the development of the scale-invariant 
MERA with a boundary described in Evenbly et al. (2010). The concept of a boundary 
scale- invariant MERA was proposed by G. Vidal, with attempts at implementation by 
R. N. C. Pfeifer and G. Evenbly, and theoretical support from V. Pico, S. Iblisdir, L. Tagli- 
acozzo, and I. McCulloch. The author's implementation attempted to optimise both the 
bulk and the boundary of the MERA simultaneously, and was not overly successful. It was 
subsequently laid aside in favour of the implementation by G. Evenbly described in the 
above reference. The interested reader is directed to this paper for further information. 
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Chapter 3 

Abelian Symmetries of Spin Systems 



In Chapter 2, I explained how it is possible to exploit the scale invariance of a quantum 
critical system on a lattice to construct an efficient Ansatz for the description of the low 
energy subspace, and how this could be used to extract the conformal data describing 
the behaviour of this system in the infra-red (large length-scale) limit. Although other 
numerical techniques and different Ansatze have previously been employed to study such 
systems, the approach developed in Chapter 2 was unique in the way in which it reflects 
the underlying scaling symmetry of the system. This resulted in an Ansatz which nat- 
urally reproduced the polynomially decaying correlators of a quantum critical system, 
and provided easy access to the data of the associated conformal field theory, as well 
as providing a numerical description of an infinite physical system which is remarkably 
compact. 

But scale invariance is not the only symmetry exhibited by quantum lattice models. 
Frequently a Hamiltonian will exhibit additional internal global symmetries which may 
be described by a group. For example, the Hamiltonian of the Ising model. 




(3.1) 



s 
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is invariant under a rotation of tt radians about the z axis, 
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giving it a Z2 symmetry. The XX model, 

Hxx = -JX: (ai^)^r^) + -hY,^i^\ (3.6) 

s=l s=l 

similarly exhibits a U(l) symmetry corresponding to invariance under rotation by any 
angle about the z axis, and the Heisenberg model, 

^Heisenberg = -J^ i^-A'^ ^i'^'^ + MI'^ + jA'^^i'^'^ + ^^i'^) ' (3-7) 

s 

has a U(l) symmetry about the z axis for = Jy (known as the Heisenberg XXZ 
model), and an SU(2) symmetry if Jx = Jy = Jz and h = (the Heisenberg XXX 
model). The natural question to ask was whether global symmetries such as these can 
also be exploited, either to facilitate the study of a particular symmetry sector of the 
model, or for computational advantage. (As an example of the former, in the Heisenberg 
XXX model, a U(l) subgroup of the SU(2) symmetry may be identified with particle 
number, and it may be desireable to study a lattice with a particular number of particles 
present. This may be achieved approximately by using a chemical potential, but the 
resulting model is still subject to fluctuations in particle number. It would be useful to 
be able to fix exactly either the particle number density or the total particle number on 
the lattice, and we shall see that this can indeed be done in Sec. 3.2.) 



In this Chapter I will therefore discuss the exploitation of global symmetries of the Hamil- 
tonian of a quantum lattice model, with particular attention to Abelian symmetries. This 
Chapter is divided into three parts: 
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1. A concise summary of how these internal symmetries manifest in tensor network 
Ansatze (Sec. 3.1). 

2. An example, being a self-contained presentation of the exploitation of U(l) symme- 
try in the MERA (Sec. 3.2). 

3. A discussion of the practicalities of how these symmetries may be efficiently imple- 
mented for an Abelian symmetry group, and how this is modified in the presence of 
fermionic exchange statistics (Sec. 3.5). 

The focus of Sec. 3.1 will be predominantly on development of the formalisms and tech- 
niques respectively whereby the global internal symmetries of a Hamiltonian may be 
exploited. Section 3.2 puts this material into practice, presenting in considerable detail 
how the formalisms of the preceding Sections are applied in the construction of a U(l)- 
invariant MERA, culminating in demonstrations of both the ability to select out any sym- 
metry sectors of a system which may be desired, and a substantial (approximately eight- 
to tenfold) decrease in computational cost when compared with the standard MERA, 
confirming that it is indeed both feasible and useful to exploit these symmetries in the 
MERA. Section 3.5 provides an additional level of implementation detail not present in 
Sec. 3.2, and discusses the extension of the approach presented here to systems of fermions, 
which have non-trivial exchange statistics in addition to a Z2 parity symmetry. 



Section 3.1 of this Chapter has previously been published as Singh, Pfeifer, and Vidal, 
Physical Review A, 82, 050301, 2010, © (2010) by the American Physical Society. 

Sections 3.2-3.4 of this Chapter have previously been published as Singh, Pfeifer, and 
Vidal, Physical Review B, 83, 115125, 2011, © (2011) by the American Physical Society. 
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3.1 Tensor Network Decompositions in the Presence 
of a Global Symmetry 

Tensor network decompositions offer an efficient description of certain many-body states 
of a lattice system and are tlie basis of a wealtli of numerical simulation algorithms. In 
this Section I discuss how to incorporate a global symmetry, given by a compact, com- 
pletely reducible group Q, into tensor network decompositions and algorithms. This is 
achieved by considering tensors that are invariant under the action of the group Q. Each 
symmetric tensor decomposes into two types of tensors: degeneracy tensors, containing 
all the degrees of freedom, and structural tensors, which depend only on the symmetry 
group. In numerical calculations, the use of symmetric tensors ensures the preservation 
of the symmetry, allows selection of a specific symmetry sector, and significantly reduces 
computational costs. On the other hand, the resulting tensor network may also be inter- 
preted as a superposition of exponentially many spin networks. Spin networks are used 
extensively in loop quantum gravity, where they represent states of quantum geometry. 
This work highlights their importance also in the context of tensor network algorithms, 
thus setting the stage for cross-fertilization between these two areas of research. 

3.1.1 Introduction 

Locality and symmetry are pivotal concepts in the formulation of physical theories. In 
a quantum many-body system, locality implies that the dynamics are governed by a 
Hamiltonian H that decomposes as the sum of terms involving only a small number of 
particles, and whose strength decays with the distance between the particles. In turn, a 
symmetry of the Hamiltonian H allows us to organize the kinematic space of the theory 
according to the irreducible representations of the symmetry group. 

Both symmetry and locality can be exploited to obtain a more compact description of 
many-body states and to reduce computational costs in numerical simulations. In the 
case of symmetries, this has long been understood. Space symmetries, such as invariance 
under translations or rotations, as well as internal symmetries, such as particle number 
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conservation or spin isotropy, divide the Hilbert space of the theory into sectors labeled 
by quantum numbers or charges. The Hamiltonian H is by definition block-diagonal in 
these sectors. If, for instance, the ground state is known to have zero momentum, it can 
be obtained by just diagonalizing the (comparatively small) zero-momentum block of H. 

In recent times, the far-reaching implications of locality for our ability to describe many- 
body systems have also started to unfold. The local character of the Hamiltonian H limits 
the amount of entanglement that low-energy states may have, and in a lattice system, 
restrictions on entanglement can be exploited to succinctly describe these states with 
a tensor network decomposition. Examples of tensor network decompositions include 
Matrix Product States (Fannes et al., 1992b; Ostlund and Rommer, 1995), Projected 
Entangled-Pair States (Sierra and Martm-Delgado, 1998; Verstraete and Cirac, 2004), 
and the MERA (Vidal, 2007a, 2008). Importantly, in a lattice made of sites, where 
the Hilbert space dimension grows exponentially with N , tensor network decompositions 
often offer an efficient description (with costs that scale roughly as A^). This allows for 
scalable simulations of quantum lattice systems, even in cases that are beyond the reach 
of standard Monte Carlo sampling techniques. As an example, the MERA has been 
recently used to investigate ground states of frustrated antiferromagnets (Evenbly and 
Vidal, 2010a). 

In this Section we investigate how to incorporate a global symmetry into a tensor net- 
work, so as to be able to simultaneously exploit both the locality and the symmetries of 
physical Hamiltonians to describe many-body states. Specifically, in order to represent 
a symmetric state that has a limited amount of entanglement, we use a tensor network 
made of symmetric tensors. This leads to an approximate, efficient decomposition that 
preserves the symmetry exactly. Moreover, a more compressed description is obtained 
by breaking each symmetric tensor into several degeneracy tensors (containing all the 
degrees of freedom of the original tensor) and structural tensors (completely fixed by 
the symmetry). This decomposition leads to a substantial reduction in computational 
costs and reveals a connection between tensor network algorithms and the formalism of 
spin networks (Penrose, 1971b) used in loop quantum gravity (Rovelli, 2008; Rovelli and 
Smolin, 1995). 
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In the case of an MPS, global symmetries have already been studied by many authors 
in the context of both one- dimensional quantum systems and two-dimensional classical 
systems (see e.g. Fannes et al. 1992b; McCuUoch and Gulacsi 2002; McCuUoch 2007; 
Ostlund and Rommer 1995; Perez-Garcia et al. 2008; Sanz et al. 2009; Sierra and Nishino 
1997; Singh et al. 2010b; White 1992). An MPS is a trivalent tensor network (i.e., each 
tensor has at most three indices) and symmetries are comparatively easy to characterize. 
The present analysis applies to the more challenging case of a generic tensor network 
decomposition (where tensors typically have more than three indices). 

3.1.2 Symmetric Decomposition of a Tensor Network 

We consider a lattice L made of sites, where each site is described by a complex vector 
space V of finite dimension d. A pure state |\E') G V®^ of the lattice can be expanded as 

d 

|^)= ^ (^)n*2...i^Kl^2...^iv), (3.8) 

il,i2vi«iV=l 

where l^^) denotes a basis of V for site s G £. For our purposes, a tensor network 
decomposition for |^) consists of a set of tensors T'^'"^ and a network pattern or graph 
characterized by a set of vertices and a set of directed edges. Each tensor T'-'') sits at a 
vertex v of the graph, and is connected with neighboring tensors by hond indices according 
to the edges of the graph. The graph also contains open edges, corresponding to the N 
physical indices 21,^2, • • • lin- The d^ coefficients {^)i^i2...iN expressed as [Fig. 3.1(i)] 




namely as the tensor product of the tensors T'^'"^ on all the vertices f , where the tensor 
trace tTr contracts all bond indices, so that only the physical indices zi, 12, . . .iN remain 
on the r.h.s. of Eq. (3.9). 

We also introduce a compact, completely reducible group Q. This includes finite groups 
as well as Lie groups such as 0(n), SO(n), U(n), and SU(n). Let U : Q ^ L(y) be a 
unitary matrix representation of Q on the space V of one site, so that for each g & Q, 
f/g : V — 7- V denotes a unitary matrix and Ug-^g^ = Ug^Ug^. Here we are interested in states 
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that are invariant under transformations of the form [7®^/ 

(f/,)«^|vl>) = |v]>), \Jgeg. (3.10) 
The space V of one site decomposes as the direct sum of irreducible representations (irreps) 

V ^ = (D'* ® V^) , (3.11) 

a a 

where denotes the irrep labeled with charge a and da is the number of times appears 
in V. We denote by a = the charge corresponding to the trivial irrep, so that V'^ = C 
and Ug = 1. In Eq. (3.11) we have also rewritten the same decomposition in terms of a 
da- dimensional degeneracy space D". We choose a local basis \i) = \a, rria) in V, where 
aa labels states within the degeneracy space (i.e. = 1, . . . , da) and labels states 
within irrep V". In this basis, Ug reads 

f/^ = 0(F®f/»). (3.12) 

a 

Recall that an operator M : V — )■ V that commutes with the group, [M, Ug] = for all 
g ^ g, decomposes as (Cornwcll, 1997) 

M = (^M" (g) r) (3.13) 

a 

(Schur's lemma). 

Our goal is to characterize a tensor network made of symmetric tensors, namely, tensors 
that are invariant under the simultaneous action of Q on all their indices. A symmetric 
tensor T with for example, two outgoing indices i and j and one incoming index k, fulfills 
[Fig. 3.1(ii)] 

j2^UgUiVg)j,,{n,,iW^g),k' = (T).,,V,,, e g, (3.14) 

ijk 

where U, V, and W denote unitary matrix representations of g. Clearly, this choice 
guarantees that Eq. (3.10) is satisfied [Fig. 3.1(iii)]. Standard group representation theory 

set of states \'^t) that transform covariantly, {Ug)'^^\'i>t) = X]t'(^s)t*'l^t')i where is a unitary 
representation of can be represented by an invariant pure state |$) oc |^'t)|t) of lattice C and 
one additional site on which the group acts with Wj. The same is true for an invariant mixed state 

p « Et with {u^r^piuir^. 
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Figure 3.1 - (i) Four-site state ^ expressed in terms of a tensor network made of three 
tensors connected according to a directed graph, (ii) Invariance of tensor T in Eq. (3.14). 
(iii) Invariance of a tensor network of symmetric tensors, Eq. (3.10). 



results (Cornwcll, 1997) imply that each symmetric tensor can be further decomposed in 
such a way that the degrees of freedom that are not fixed by the symmetry can be isolated 
(Fig. 3.2). Next we discuss the cases of tensors with a small number of indices. Recall that 
an index i of a tensor is associated with a vector space that decomposes as in Eq. (3.11); 
therefore we can write i = (a, m^), j = (6, f3h, rih), k = (c, 7c, Oc), and so on. 

One leg. — A tensor T with only one index i is invariant only if Q acts on it trivially, so 
the only relevant irrep is a = 0, and index i = labels states within the degeneracy 
space V°. 

Two legs. — Schur's lemma (Cornwell, 1997) establishes that a symmetric tensor T with 
one outgoing index i and one incoming index j decomposes as [cf. Eq. (3.13)] 

{T)ij = {P°''')aal3biQ"''')mant,, Q"^ = ^ab^man^- (3.15) 

Thus, for fixed values of the charges a and 6, {T)ij breaks into a degeneracy tensor P^^ 
(where only a = 6 is relevant) and another tensor Q"-^ . P""^ contains all the degrees of 
freedom of T that are not fixed by the symmetry, whereas Q""^ is completely determined 
by Q. Another combination of outgoing and incoming indices, for example two incoming 
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Figure 3.2 - Decomposition of tensors with one to four indices. The sums in (iv) run over 
the intermediate indices {e,ee,qe) and (/, C/i^/) Eqs. (3.18)-(3.19). 

indices, leads to a different form for tensor Q"*. 

Three legs. — The tensor product of two irreps with charges a and b can be decomposed 
as the direct sum of irreps, 

V" (g) ^ iV^^feV^ (3. 16) 

c 

where A^"^^ denotes the number of copies of that appear in the tensor product. For 
notational simplicity, from now on we assume that Q is multiplicity-free,^ that is, A*"^^ < 1, 
and denote by {Q"'^'^) manioc the change of basis between the product basis \a,ma) ® \b,nh) 
and the coupled basis |c, Oc). The Wigner-Eckart theorem states that a symmetric tensor 
T with, for example, two outgoing indices i,j and one incoming index k, then decomposes 
as 

(3.17) 

As before, for fixed values of the charges a, b, and c, {T)ijk factorizes into degeneracy 
tensors P"*"^ with all the degrees of freedom and structural tensors Q"'^'^ (the Clebsch- 
Gordan coefficients) completely determined by the group Q. An analogous decomposition 
with different Q°'^'^ holds for other combinations of incoming and outgoing indices. 

Four legs. — The tensor product of three irreps V" (g) ® may contain several copies of 
an irrep V^. Let e be the charge that results from fusing a and b, V* (g) = 0^ A'^^V^. 

^In non-multiplicity-free groups, such as SU(3), where N^^^ might be larger than 1, the coupled basis 
|c, Oc, /i) and tensor S'^'"^ must include an extra index /x = 1, . . . , A^^j,. See, for example. Chapter 4 of this 
Thesis, and Pfcifcr ct al. (2010). 
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We can use the values of e for which N^j^N^^ ^ (i.e., such that a and h fuse to e, and 
e and c fuse to d) to label the different copies of that appear in V" (g) (g) V^. Let 
{Qf"^'^)manbOcPdqe denote the change of basis between the product basis \ama) ® \bnb) |cOc) 
and the coupled basis \dpd]e) obtained by fusing to the intermediate basis \eqe) G V^. 
Then a symmetric tensor T with three outgoing indices i, j, and k and one incoming 
index / = {d,Sd,Pd) decomposes as'^ 

{Qe) (3.18) 

e 

where the sum is over all relevant values of the intermediate charge e. Alternatively, T 
can be decomposed as 

{T)ijkl = 'y^S^f^'^'^)aal3fycSd{Qf"^'^)mantOcPd^ (3.19) 
/ 

where {Qf"^'^)mantOcPd denotes the change of basis to another coupled basis \dpd; f) of 
obtained by fusing first b and c into /, and then a and / into d, involving a different inter- 
mediate charge index /. The two coupled bases are related by a unitary transformation 
given by the 6-index tensor F [related to the 6-j symbols for e.g. Q = SU(2); see Eq. (5.7) 
of Chapter 5] such that 

Qabcd ^ J2{FfyjQf^'^. (3.20) 

e 

Since Eqs. (3.18) and (3.19) represent the same tensor T, the degeneracy tensors P and 
P are related by 

pabcd ^ y^^^pabc*y^pabcd_ {S.21) 

e 

More generally, a symmetric tensor T with t indices is = (a^, a^^, "niaj, where s = 1, . . . ,t, 
decomposes as 

{T)iii2...it = ^^(-Pei^.ej!)aai---aat (Qei...e*/)mai---ma( ; (3.22) 

^Note that Eq. (3.18) differs from Eq. (11) of Singh et al. (2010a). This is due to an error in the 
pubhshed paper, where this equation was mistakenly given as 

{T)ijkl — y ] {Pe )aapblcSdta{Qe )marn,OaPdq^J (H) 

with exphcit degeneracy indices eg and qe associated with the intermediate charge index e. For a given 
value of Ee, {Pe^'^'^)aal3biaSdt„ nou-zcro for precisely one set of values {aa, /?(,, 7c, ^d}, and similarly for 
1e: iQe'"^'^)mantOaPdqa^ ^.nd {TTia, nf,, Oc, Pd} ■ Consequently, to obtain the most efficient representation of 
(T)ijki we would always evaluate the sum over and qe in the above expression, reducing it to Eq. (3.18). 
Similar corrections have been made to Eqs. (3.19) and (3.22). 
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Figure 3.3 - A tensor network for a symmetric state l^*) G v®^ of lattice £ (Fig. 3.1) 
is expressed as a linear superposition of spin networks. The sum runs over the intermediate 
indices that carry charges e and / (shown explicitly) as well as all indices shared by two tensors. 



where the sum is over the intermediate charges Ck, k = 1, . . . ,t' . The degeneracy tensors 
Pel"'e* contain all the degrees of freedom of T, whereas the structural tensors Qei.'.'.'e*/ 
are completely determined by the symmetry. Here ei,e2, ■ ■ ■ ,et' are intermediate charges 
that decorate the inner branches of a trivalent tree used to label a basis in the space of 
intertwining operators between the tensor products of incoming and outgoing irreps. A 
different choice of tree will produce different sets of tensors P and Q, related to P and Q 
by F- moves.'' 

We can now investigate how the tensor network decomposes if we write each of its tensors 
T in the (P, Q) form of Eq. (3.22) (see Fig. 3.3). For any fixed value of all the charges, the 
whole tensor network factorizes into two terms. The first one is a network of degeneracy 
tensors. The second one is a directed graph with edges labeled by irreps of Q and vertices 
labeled by intertwining operators. This is nothing other than a spin network (Penrose, 
1971b), a well-known object in mathematical physics and, especially, in loop quantum 
gravity (Rovelli, 2008; Rovelli and Smolin, 1995), where it is used to describe states of 
quantum geometry. Accordingly, a symmetric tensor network for the state |\1') G V®^ of a 
lattice C of N sites can be regarded as a linear superposition of spin networks with N open 
edges. The number of spin networks in the linear superposition grows exponentially with 
the size of the tensor network. The expansion coefficients are given by the degeneracy 
tensors. 

"^When Q is an Abelian group, such as U(l), the tensor product V of two hreps only gives rise to 
one irrep V^, so that no intermediate charges ei, 62, ... , ec need to be specified in Eq. (3.22), simphfying 
significantly the decomposition of symmetric tensors. 
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= AA = A" 

Figure 3.4 - Product of two symmetric tensors. Only the intermediate charges d, e, and / are 
explicitly shown. Additional sums apply to all indices shared by two tensors. The computation 
involves evaluating spin networks. 




3.1.3 Applications of Symmetric Tensor Networks 

Computationally, the present characterization of a symmetric tensor network is of interest 
for several reasons. First, it allows us to describe a state with specific quantum 

numbers, which are preserved exactly during approximate numerical simulations. Let us 
consider as an example the group U(l), with charge n corresponding to particle number 
{n = 0, ±1, ±2, . . .), and the group SU(2), with charge j corresponding to the spin (j = 
0, 1/2, 1, 3/2, . . .). The symmetric tensor network can be used to describe a state with, for 
example, zero particles (n = 0) and zero spin (j = 0), respectively — or, more generally, 
covariant states with any value of n and j.^ 

Second, the (P, Q)-decomposition (3.22) concentrates all the degrees of freedom of a sym- 
metric tensor T in the degeneracy tensors P, producing a more compact description. For 
instance, for the U(l) and SU(2) groups, an approximation of the ground state of the an- 
tiferromagnetic Heisenberg spin-| chain with a MERA of bond dimension x = 21 requires 
five and thirty-five times less parameters than with nonsymmetric tensors, respectively 
(Singh and Vidal, in preparation; Singh et al., 2011). 

In addition, the (P, (5)-decomposition (3.22) lowers the cost of simulations significantly. 
Consider the multiplication of two tensors (Fig. 3.4) which is central to most tensor 
network algorithms. Cost reductions come from two fronts: 

^See footnote 1, above. 
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1. Block-sparse matrices. — The most costly step in multiplying two tensors T' and T" 
consists of multiplying two matrices M' and M" obtained from T' and T". These 
matrices are of the form of Eq. (3.13), and therefore their multiplication can be done 
blockwise: 

M = M'M" = [(M"'M"") r . (3.23) 

a 

2. Pre- computation. — Given a (P, (5)-decomposition of tensor T, another (P, (5)-decomposition 
(as required, e.g., to obtain the matrices M' and M" above) involves a linear map 

T: 

P = r(P). (3.24) 

This map F, of which Eq. (3.21) is an example, is completely determined by the 
symmetry. In those tensor network algorithms that proceed by repeating a sequence 
of manipulations, map F can be computed once and stored in memory for repeated 
usage. 

More detailed explanations of algorithmic details, as well as practical examples of the 
gains obtained using invariant tensors, are presented in Singh et al. (2011) (also Sec. 3.2 
of this Thesis) and Singh and Vidal (in preparation) for the groups U(l) and SU(2), 
respectively. Evenbly and Vidal (2010a) exploited the U(l) symmetry in a 2D MERA 
calculation that involved tensors with up to twelve indices. 

Finally, the connection between symmetric tensor networks and spin networks allows us 
to import into the context of tensor network algorithms techniques developed to evaluate 
spin networks in loop quantum gravity. Such techniques can be used, for example, to 
compute the linear map F of Eq. (3.24). Conversely, tensor network algorithms may 
also prove useful in loop quantum gravity, since they allow (for example) the efficient 
manipulation of superpositions of an exponentially large number of spin networks. 



3.2 Example: U(l)-Symmetric MERA 

Tensor network decompositions offer an efficient description of certain many-body states 
of a lattice system and are the basis of a wealth of numerical simulation algorithms. In 
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Sec. 3.1 I discussed how to incorporate a global internal symmetry, given by a compact, 
completely reducible group Q, into tensor network decompositions and algorithms. Here 
I specialize to the case of Abelian groups and, for concreteness, to a U(l) symmetry, 
associated, e.g., with particle number conservation. I will consider tensor networks made 
of tensors that are invariant (or covariant) under the symmetry, and explain how to de- 
compose and manipulate such tensors in order to exploit their symmetry. In numerical 
calculations, the use of U(l)-symmetric tensors allows selection of a specific number of 
particles, ensures the exact preservation of particle number, and significantly reduces com- 
putational costs. I illustrate all these points in the context of the multi-scale entanglement 
renormalization Ansatz. 

3.2.1 Introduction 

Tensor networks are becoming increasingly popular as a tool to represent wave-functions of 
quantum many-body systems. Their success is based on the ability to efficiently describe 
the ground state of a broad class of local Hamiltonians on the lattice. Tensor network 
states are used both as a variational Ansatz to numerically approximate ground states 
and as a theoretical framework to characterize and classify quantum phases of matter. 

Examples of tensor network states for one dimensional systems include the Matrix Prod- 
uct State or MPS*^, which results naturally from both Wilson's numerical renormalization 
group (Wilson, 1975) and White's Density Matrix Renormalization Group (DMRG),' and 
is also used as a basis for simulation of time evolution, e.g. with the time evolving block 
decimation (TEBD) algorithm (Vidal, 2007b, 2003, 2004) and variations thereof, often 
collectively referred to as time-dependent DMRG;*^ the Tree Tensor Network (TTN) (Shi 
et al., 2006), which follows from coarse-graining schemes where the spins are blocked hi- 
erarchically; and the Multi-scale Entanglement Renormalization Ansatz (MERA),'^ which 
results from a renormalization group procedure known as entanglement renormalization 
(Vidal, 2007a, 2010). For two dimensional lattices there are generalizations of these three 

•^Fannes et al. (1992a); Ostlund and Rommer (1995); Perez-Garcia et al. (2007) 
^McCuUoch (2008); SchoUwock (2005a, 2011); White (1992, 1993) 

^Daley et al. (2004); Schollwock (2005b); Vidal (2007b, 2003, 2004); White and Feiguin (2004) 
^Evenbly and Vidal (2009a); Giovannetti ct al. (2008); Pfeifer et al. (2009); Vidal (2007a, 2008, 2010) 
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tensor network states, namely projected entangled pair states (PEPS)/" 2D TTN (Murg 
et al., 2010; Tagliacozzo ct al., 2009), and 2D MERA,^^ respectively. As variational 
Ansatze, PEPS and 2D MERA are particularly interesting since they can be used to ad- 
dress large two-dimensional lattices, including systems of frustrated spins (Evenbly and 
Vidal, 2010a; Murg et al., 2009) and interacting fermions,^^ where Monte Carlo techniques 
fail due to the sign problem. 

A many-body Hamiltonian H may be invariant under transformations that form a group of 
symmetries (Cornwell, 1997). The symmetry group divides the Hilbert space of the theory 
into symmetry sectors labeled by quantum numbers or conserved charges. On a lattice 
one can distinguish between space symmetries, which correspond to some permutation of 
the sites of the lattice, and internal symmetries, which act on the vector space of each site. 
An example of space symmetry is invariance under translations by some unit cell, which 
leads to conservation of momentum. An example of internal symmetry is SU(2) invariance, 
e.g. spin isotropy in a quantum spin model. An internal symmetry can in turn be global, 
if it transforms the space of each of the lattice sites according to the same transformation 
(e.g. a spin-independent rotation); or local, if each lattice site is transformed according to 
a different transformation (e.g. a spin- dependent rotation), as it is for gauge symmetric 
models. A global internal SU(2) symmetry gives rise to conservation of total spin. By 
targetting a specific symmetry sector during a calculation, computational costs can often 
be significantly reduced while explicitly preserving the symmetry. It is therefore not 
surprising that symmetries play an important role in numerical approaches. 

In Sec. 3.1 I described a formalism for incorporating global internal symmetries into a 
generic tensor network algorithm. Both Abelian and non-Abelian symmetries were con- 
sidered. The purpose of Sec. 3.2 is to address, at a pedagogical level, the implementation 
of Abelian symmetries into tensor networks. We will also discuss several more practical 
aspects of the exploitation of Abelian symmetries not covered in Sec. 3.1. For concreteness 

i°Gu et al. (2008); Jiang et al. (2008); Jordan et al. (2008); Murg et al. (2007, 2009); Nishino and 
Okunishi (1998); Nishio ct al. (2004); Sierra and Martin-Delgado (1998); Verstraete and Cirac (2004); 
Xie et al. (2009) 

"Aguado and Vidal (2008); Cincio et al. (2008); Evenbly and Vidal (2009b, 2010a,b,c); Konig et al. 
(2009) 

i^Barthel et al. (2009); Corboz and Vidal (2009); Corboz ct al. (2010a,b); Gu et al. (2010); Kraus et al. 
(2010); Li et al. (2010); Pineda et al. (2010); Pizorn and Verstraete (2010); Shi et al. (2009) 
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this Section concentrates on U(l) symmetry, but extending these results to any Abehan 
group is straightforward. A similar analysis of non-Abelian groups will be considered in 
Singh and Vidal (in preparation), as well as being discussed briefly in Chapter 5 of this 
Thesis. 

In tensor network approaches, the exploitation of global internal symmetries has a long 
history, especially in the context of MPS. Both Abelian and non-Abelian symmetries 
have been thoroughly incorporated into DMRG code and have been exploited to obtain 
computational gains. ^'^ Symmetries have also been used in more recent proposals to 
simulate time evolution with MPS.^^ 

When considering symmetries, it is important to notice that an MPS is a trivalent tensor 
network. That is, in an MPS each tensor has at most three indices. The Clebsch-Gordan 
coefficients, or coupling coefficients, of a symmetry group are also trivalent (Cornwell, 
1997), and this makes incorporating the symmetry into an MPS by considering symmetric 
tensors particularly simple. In contrast, tensor network states with a more elaborate 
network of tensors, such as MERA or PEPS, consist of tensors having a larger number of 
indices. In this case a more general formalism is required in order to exploit the symmetry. 
As explained in Sec. 3.1, a generic symmetric tensor can be decomposed into a degeneracy 
part, which contains all degrees of freedom not determined by symmetry, and a structural 
part, which is completely determined by symmetry and can be further decomposed as a 
trivalent network of Clebsch-Gordan coefficients. 

The use of symmetric tensors in more complex tensor networks has also been discussed 
in Pcrcz-Garcia et al. (2010) and Zhao et al. (2010). In particular, Perez-Garcia et al. 
(2010) has shown that under convenient conditions (injectivity), a PEPS that represents a 
symmetric state can be represented with symmetric tensors, generalizing similar results for 
MPS obtained in Pcrcz-Garcia et al. (2008). Notice that these studies are not concerned 
with how to decompose symmetric tensors so as to computationally protect or exploit 

i^Bergkvist et al. (2006); McCuUoch and Gulacsi (2002); McCuUoch (2007); Ostlund and Rommer 
(1995); Perez-Garcfa et al. (2008); Pittcl and Sandulescu (2006); Ramasesha et al. (1996); Sanz ct al. 
(2009); Schollwock (2005b); Sierra and Nishino (1997); Tatsuaki (2000); White (1992) 

"Cai et al. (2010); Daley ct al. (2004, 2005); Danshita et al. (2007); Mishmash et al. (2009); Muth 
et al. (2010); Schollwock (2005b); Singh et al. (2010b); Vidal (2007b, 2004); White and Fciguin (2004) 
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the symmetry. On the other hand, exploitation of U(l) symmetry for computational gain 
in the context of PEPS was reported in Zhao et al. (2010), although no implementation 
details were provided. Finally, several aspects of local internal symmetries in tensor 
network algorithms have been addressed in Schuch ct al. (2010), Swingle and Wen (2010), 
Chen et al. (2010b), and Tagliacozzo and Vidal (2010). 

The discussion of the U(l)-symmetric MERA is organized into Sections as follows: 

Section 3.2.2 contains a review of the tensor network formalism and introduces the nomen- 
clature and diagrammatical representation of tensors used in the rest of the Chapter. It 
also describes a set V of primitives for manipulating tensor networks, consisting of ma- 
nipulations that involve a single tensor (permutation, fusion and splitting of the indices 
of a tensor) and matrix operations (multiplication and factorization). 

Section 3.2.3 reviews basic notions of representation theory of the Abelian group U(l). 
The action of the group is analysed first on a single vector space, where U(l)-symmetric 
states and U(l)-invariant operators are decomposed in a compact, canonical manner. This 
canonical form allows us to identify the degrees of freedom which are not constrained by 
the symmetry. The action of the group is then also analysed on the tensor product of two 
vector spaces and, finally, on the tensor product of a finite number of vector spaces. 

Section 3.2.4 explains how to incorporate the U(l) symmetry into a generic tensor network 
algorithm, by considering U(l)-invariant tensors in a canonical form, and by adapting the 
set V of primitives for manipulating tensor networks. These include the multiplication of 
two U(l)-invariant matrices in their canonical form, which is at the core of the computa- 
tional savings obtained by exploiting the symmetry in tensor network algorithms. 

Section 3.2.5 illustrates the practical exploitation of the U(l) symmetry in a tensor net- 
work algorithm by presenting MERA calculations of the ground state and low energy 
states of two quantum spin chain models. 

The canonical form offers a more compact description of U(l)-invariant tensors, and leads 
to faster matrix multiplications and factorizations. However, there is also an additional 
cost associated with maintaining an invariant tensor in its canonical form while reshaping 
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(fusing and/or splitting) its indices. In some situations, this cost may offset the benefits 
of using the canonical form. In Sec. 3.4 we discuss a scheme to lower this additional cost 
in tensor network algorithms that are based on iterating a repeated sequence of trans- 
formations. This is achieved by identifying, in the manipulation of a tensor, operations 
which only depend on the symmetry. Such operations can be precomputed once at the 
beginning of a simulation. Their result, stored in memory, can be re-used at each iteration 
of the simulation. Section 3.4 describes two such specific precomputation schemes. 

3.2.2 Review: Tensor Network Formalism 

In this Section we review background material concerning the formalism of tensor net- 
works, without reference to symmetry. We introduce basic definitions and concepts, as 
well as the nomenclature and graphical representation for tensors, tensor networks, and 
their manipulations, that will be used in Sec. 3.2. 

3.2.2.1 Tensors 

A tensor T is a multidimensional array of complex numbers Ti^i^...i^. G C. The rank of 
tensor T is the number k of indices. For instance, a rank-0 tensor {k = 0) is a complex 
number. Similarly, rank-1 {k = 1) and rank-2 [k = 2) tensors represent vectors and 
matrices, respectively. The size of an index i, denoted is the number of values that 
the index takes, i G {1, 2, . . . , The size of a tensor T, denoted |T|, is the number of 
complex numbers it contains, namely |T| = x |i2| x . . . x In this discussion of 
the U(l)-symmetric MERA, we will use the hat (") to indicate that an object is a tensor. 
Vectors are included in this convention, writing their components as, e.g., although 
for simplicity we will omit the hat when a vector is written in bra or ket form, e.g. |\E'). 

It is convenient to use a graphical representation of tensors, as introduced in Fig. 3.5, 
where a tensor T is depicted as a circle (more generally some shape, e.g. a square) and 
each of its indices is represented by a line emerging from it. In order to specify which index 
corresponds to which emerging line, we follow the prescription that the lines corresponding 
to indices {ii,i2, ■ ■ ■ ,ik} emerge in counterclockwise order. Unless stated otherwise, the 
first index will correspond to the line emerging at nine o'clock (or the first line encountered 
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(0 




h 




I 



h 



Figure 3.5 - (i) Graphical representation of a tensor T of rank k and components Tij^i2...ik- 
The tensor is represented by a shape (circle) with k emerging lines corresponding to the k 
indices ii,i2, ■ ■ ■ , ik- Notice that the indices emerge in counterclockwise order, (ii) Graphical 
representation of tensors with rank A; = 0, 1, and 2, corresponding to a complex number 
c G C, a vector \v) G Cl*l, and a matrix M G Cl^il'^l^^l, respectively. 

while proceeding counterclockwise from nine o'clock). 

Two elementary ways in which a tensor T can be transformed are by permuting and 
reshaping its indices. A permutation of indices corresponds to creating a new tensor T' 
from T by simply changing the order in which the indices appear, e.g. 



On the other hand, a tensor T can be reshaped into a new tensor T' by "fusing" and/or 
"splitting" some of its indices. For instance, in 



tensor T' is obtained from tensor T by fusing indices 6g{1,...,|6|} and cG{1,...,|c|} 
together into a single index d of size \d\ = \b\ ■ \c\ that runs over all pairs of values of b 
and c, i.e. d G {(1, 1), (1, 2), . . . , |c| — 1), |c|)}, whereas in 



(3.25) 



d = b X C, 



(3.26) 



Tabc (T ^adi 



d = b X c, 



(3.27) 



tensor T is recovered from T' by splitting index d of T' back into indices b and c. The 
permutation and reshaping of the indices of a tensor have a straightforward graphical 
representation; see Fig. 3.6. 
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^ b 

Figure 3.6 - Transformations of a tensor: (i) Permutation of indices b and c. (ii) Fusion of 
indices b and c into d = b x c; splitting of index d = b x c into 6 and c. 



T R S 




Figure 3.7 - (i) Graphical representation of the multiplication of two matrices R and S to 
give a new matrix T (3.28) (ii) Graphical representation of an example contraction of two 
tensors R and S into a new tensor T (3.29). 



3.2.2.2 Multiplication of Two Tensors 

Given two matrices R and S with components Rat and She, we can multiply them together 
to obtain a new matrix T, T = R ■ S , with components 

Tac = ""^^ RabSbc, (3.28) 
b 

by summing over or contracting index b. The multiplication of matrices R and S is rep- 
resented graphically by connecting together the emerging lines of R and S corresponding 
to the contracted index, as shown in Fig. 3.7(i). 

Matrix multiplication can be generalized to tensors. For instance, given tensors R and S 
with components Rabcd and Scfbh, we can define a tensor T with components Thafd given 
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Figure 3.8 - Graphical representations of the five elementary steps 1-5 into which one can 
deconnpose the contraction of the tensors of Eq. (3.29). 



by 

Thafd = 'y^^RgbcdScfbh- (3.29) 
be 

Again the multiplication of two tensors can be graphically represented by connecting 
together the lines corresponding to indices that are being contracted [indices b and c in 
Eq. (3.29)]; see Fig. 3.7(ii). 

The multiplication of two tensors can be broken down into a sequence of elementary steps 
by transforming the tensors into matrices, multiplying the matrices together, and then 
transforming the resulting matrix back into a tensor. These steps are now described for 
the contraction given in Eq. (3.29). They are illustrated in Fig. 3.8. 



1. Permute the indices of tensor R in such a way that the indices to be contracted, b 
and c, appear in the last positions and in a given order, e.g. be; similarly, permute 
the indices of S so that the indices to be contracted, again b and c, appear in the 
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first positions and in the same order be: 

ad be Rabcd 

{S')bc fh = Scfbh (3.30) 

2. Reshape tensor R' into a matrix R" by fusing into a single index u all the indices 
that are not going to be contracted, u = axd, and into a single index y all indices to 
be contracted, y = b x c. Similarly, reshape tensor S' into a matrix S" with indices 
y = b X c and w = f x h, 

adbc 

{S")yw = iS')bcfh- (3.31) 

3. Multiply matrices R" and S" to obtain a matrix T", with components 

(T )uw = ^ ^ (-^ )uy {S )yw (3.32) 

y 

4. Reshape matrix T" into a tensor T' by splitting indices u = a x d and w = f x h, 

{T')adfh = (7"%«, (3.33) 

5. Permute the indices of T' into the order in which they appear in T, 

Thafd = {T')adfh- (3.34) 



Note that breaking down a multiplication of two tensors into elementary steps is not 
necessary — one can simply implement the contraction of Eq. (3.29) as a single process. 
However, it is often more convenient to compose the above elementary steps since, for 
instance, in this way one can use existing linear algebra libraries for matrix multiplica- 
tion. In addition, it can be seen that the leading computational cost in multiplying two 
large tensors is not changed when decomposing the contraction in the above steps. In 
Sec. 3.2.4.9 this subject will be discussed in more detail for U(l)-invariant tensors. 
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Figure 3.9 - (i) Factorization of a matrix T according to a singular value decomposition 
(3.35). (ii) Factorization of a rank-4 tensor T according to one of several possible singular 
value decompositions. 



3.2.2.3 Factorization of a Tensor 

A matrix T can be factorized into the product of two (or more) matrices in one of several 
canonical forms. For instance, the singular value decomposition 

Tab = ^ UacScdVdb = ^ UacScVcb (3.35) 
c,d c 

factorizes T into the product of two unitary matrices U and V, and a diagonal matrix 
5* with non-negative diagonal elements Sc = Sec known as the singular values of T; see 
Fig. 3.9(i). On the other hand, the eigenvalue or spectral decomposition of a square matrix 
T is of the form 

fab = Yl MacDcd{M-^)db = J2 MaMM~')cb, (3.36) 

c,d c 

where M is an invertible matrix whose columns encode the eigenvectors | Ac) of T, 

f|A,) = Ae|Ae), (3.37) 

M-^ is the inverse of M, and Z) is a diagonal matrix, with the eigenvalues Ac = -Dec on 
its diagonal. Other useful factorizations include the LU decomposition, the QR decom- 
position, etc. We refer to any such decomposition generically as a matrix factorization. 

A tensor f with more than two indices can be converted into a matrix in several ways, 
by specifying how to join its indices into two subsets. After specifying how tensor T is 
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Figure 3.10 - (i) Example of a tensor network M. (ii) Tensor T of which the tensor network 
Af could be a representation, (iii) Tensor T can be obtained from Af through a sequence of 
contractions of pairs of tensors. Shading indicates the two tensors to be multiplied together 
at each step. 



to be regarded as a matrix, we can factorize T according to any of the above matrix fac- 
torizations, as illustrated in Fig. 3.9(ii) for a singular value decomposition. This requires 
first permuting and reshaping the indices of T to form a matrix, then decomposing the 
latter, and finally restoring the open indices of the resulting matrices into their original 
form by undoing the reshapes and permutations. 

3.2.2.4 Tensor Networks and Their Manipulation 

A tensor network is a set of tensors whose indices are connected according to a network 
pattern, e.g. Fig. 3.10. 

Given a tensor network A/", a single tensor T can be obtained by contracting all the indices 
that connect the tensors in JV [Fig. 3.10(ii)]. Here, the indices of tensor T correspond to 
the open indices of the tensor network A/". We then say that the network A/" is a tensor 
network decomposition of T. One way to obtain T from A^ is through a sequence of 
contractions involving two tensors at a time [Fig. 3. 10 (iii)]. 
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From a tensor network decomposition M for a tensor T, another tensor network decompo- 
sition for the same tensor T can be obtained in many ways. One possibihty is to replace 
two tensors in M with the tensor resulting from contracting them together, as is done in 
each step of Fig. 3.10(iii). Another way is to replace a tensor in M with a decomposition 
of that tensor (e.g. with a singular value decomposition). In this Chapter, we will be 
concerned with manipulations of a tensor network that, as in the case of multiplying two 
tensors or decomposing a tensor, can be broken down into a sequence of operations from 
the following list: 

1. Permutation of the indices of a tensor, Eq. (3.25). 

2. Reshape of the indices of a tensor, Eqs. (3.26)-(3.27). 

3. Multiplication of two matrices, Eq. (3.28). 

4. Decomposition of a matrix [e.g. singular value decomposition (3.35) or spectral 
decomposition (3.36)]. 

These operations constitute a set V of primitive operations for tensor network manipula- 
tions (or, at least, for the type of manipulations we will be concerned with). 

In Sec. 3.2.4 we will discuss how this set V of primitive operations can be generalized to 
tensors that are symmetric under the action of the group U(l). 

3.2.2.5 Tensor Network States for Quantum Many-Body Systems 

As mentioned in Sec. 3.2.1, tensor networks are used as a means to represent the wave- 
function of certain quantum many-body systems on a lattice. Let us consider a lattice C 
made of L sites, each described by a complex vector space V of dimension d. A generic 
pure state |\E') G V®'^ of C can always be expanded as 

l^)= ^^.^2■■■^LW)\^2)■■■\^L). (3.38) 

where is = 1, . . . ,d labels a basis \is) of V for site s G £. Tensor with components 
^ni2...iL) contains complex coefficients. This is a number that grows exponentially with 
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Figure 3.11 - Examples of tensor network states for ID systems: (i) matrix product state 
(MPS), (ii) tree tensor network (TTN), (iii) multi-scale entanglement renormalization Ansatz 
(MERA). Examples of tensor network states for 2D systems: (iv) projected entangled-pair 
state PEPS, (v) 2D TTN. (2D MERA not depicted). 

the size L of the lattice. Thus, the representation of a generic pure state |^) G V®^ is 
inefficient. However, it turns out that an efficient representation of certain pure states 
can be obtained by expressing tensor \1/ in terms of a tensor network. 

Figure 3.11 shows several popular tensor network decompositions used to approximately 
describe the ground states of local Hamiltonians H of lattice models in one or two spatial 
dimensions. The open indices of each of these tensor networks correspond to the indices 
ii,Z2; ■ ■ ■ i^L of tensor ^. Notice that all the tensor networks of Fig. 3.11 contain 0{L) 
tensors. If p is the rank of the tensors in one of these tensor networks, and x is the size 
of their indices, then the tensor network depends on 0{Lx^) complex coefficients. For a 
fixed value of x this number grows linearly in L, and not exponentially. It therefore does 
indeed offer an efficient description of the pure state |^) G that it represents. Of 
course only a subset of pure states can be decomposed in this way. Such states, often 
referred to as tensor network states, are used as variational Ansatze, with the 0{Lx^) 
complex coefficients as the variational parameters. 
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Given a tensor network state, a variety of algorithms'^ are used for tasks such as: (i) com- 
putation of the expectation value (\l/|6|\l/) of a local observable 6, (ii) optimization of the 
variational parameters so as to minimize the expectation value of the energy (\E'|if|^), or 
(iii) simulation of time evolution, e.g. e~'^*|\l/). These tasks are accomplished by manip- 
ulating tensor networks. 

On most occasions, all required manipulations can be reduced to a sequence of primitive 
operations in the set V introduced in Sec. 3.2.2.4. Thus, in order to adapt tensor network 
algorithms such as those listed in footnote 15, above, to the presence of a symmetry, we 
only need to modify the set V of primitive tensor network operations. This will be done 
in Sec. 3.2.4. 

3.2.2.6 Tensors as Linear Maps 

A tensor can be used to define a linear map between vector spaces in the following way. 
First, notice that an index i can be used to label a basis of a complex vector space 
yH ^ ([AA of dimension On the other hand, given a tensor T of rank A;, we can attach 
a direction "in" or "out" to each index ii,i2, . . . ,ik- This direction divides the indices of 
T into the subset / of incoming indices and the subset O of outgoing indices. We can 
then build input and output vector spaces given by the tensor product of the spaces of 
incoming and outgoing indices, 

yti"] = (g)Vt^'', V[°"*l = 0Y^''\ (3.39) 

and use tensor T to define a linear map between V'™' and V[°"*]. For instance, if a rank-3 
tensor Tate has one incoming index c G / and two outgoing indices a,b ^ O, then it defines 
a linear map f : Vl'l ^ Vt"! ® V'^l given by 

f = J2Tabc\a)\b){c\ (3.40) 

a,b,c 

i^see e.g. Aguado and Vidal (2008); Barthel ct al. (2009); Cincio et al. (2008); Corboz and Vidal (2009); 
Corboz et al. (2010a,b); Daley et al. (2004); Evenbly and Vidal (2009a,b, 2010a,b,c); Giovannetti et al. 
(2008); Gu et al. (2008, 2010); Jiang et al. (2008); Jordan et al. (2008); Konig et al. (2009); Kraus et al. 
(2010); Li et al. (2010); MeCuUoch (2008); Murg et al. (2007, 2009, 2010); Nishino and Okunislii (1998); 
Nishio et al. (2004); Pfeifer et al. (2009); Pineda et al. (2010); SchoUwock (2005a,b, 2011); Shi et al. 
(2009, 2006); Sierra and Martm-Delgado (1998); Tagliaeozzo et al. (2009); Verstraete and Cirac (2004); 
Vidal (2007a,b, 2008, 2003, 2004, 2010); White (1992, 1993); White and Feiguin (2004); Wilson (1975); 
Xie et al. (2009). 
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Figure 3.12 - (i) Tensor T with one incoming index and two outgoing indices, denoted by 
incoming and outgoing arrows respectively (3.40). (ii) A tensor network M with directed links 
can be interpreted as a linear map between incoming and outgoing spaces (of the incoming 
and outgoing indices) obtained by composing the linear maps associated with each of the 
tensors in M. 



Graphically, we denote the direction of an index by means of an arrow; see Fig. 3.12(i). 

By decorating the lines of a tensor network Af with arrows [Fig. 3.12(ii)], this can be 
regarded as a composition of linear maps — namely, one linear map for each tensor in A/". 
While arrows might be of limited relevance in the absence of a symmetry, they will play 
an important role when we consider symmetric tensors since they specify how the group 
acts on each index of a given tensor. 



3.2.3 Review: Representation Theory of the Group U(l) 

In this Section we review basic background material concerning the representation theory 
of the group U(l). We first consider the action of U(l) on a vector space V, which 
decomposes into the direct sum of (possibly degenerate) irreducible representations. We 
then consider vectors of V that are symmetric (invariant or covariant) under the action 
of U(l), as well as linear operators that are U(l)-invariant. Then we consider the action 
of U(l) on the tensor product of two vector spaces, and its generalization to the tensor 
product of an arbitrary number of vector spaces. 
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3.2.3.1 Decomposition Into Direct Sum of Irreducible Representations 

Let V be a finite- dimensional space and let G [0, 2%) label a set of linear transformations 

W^:Y^Y, (3.41) 
that are a unitary representation of the group U(l). That is, 

WjW^ = W^Wj = 1, V e [0, 27r), (3.42) 

W^,W^, = W^,W^, = W^,+^,\,^ V ^2 e [0, 2n). (3.43) 

Then V decomposes as the direct sum of (possibly degenerate) one- dimensional irreducible 
representations (or irreps) of U(l), 

V = 0V„, (3.44) 

n 

where V„ is a subspace of dimension dn, made of (i„ copies of an irrep of U(l) with charge 
n G Z. We say that irrep n is (i„-fold degenerate and that V„ is the degeneracy space. For 
concreteness, in this Section we identify the integer charge n as labelling the number of 
particles (another frequent identification is with the z component of the spin, in which case 
semi-integer numbers may be considered). The representation of group U(l) is generated 
by the particle number operator n, 

d„ 

n = '^nPn, Pn = '^\ntn){ntn\, (3.45) 

n in=l 

where Pn is a projector onto the subspace V„ of particle number n, and the vectors |nt„), 

h\ntn) = n\ntn) , tn = l,...,dn, (3.46) 

are an orthonormal basis of V^. In terms of h, the transformations read 

= e"'^^. (3.47) 

It then follows from Eq. (3.46) that 

W^\ntn) = e-'"'^|nt„) V ^ G [0, 27r). (3.48) 
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The dual basis {{ntn\} is transformed by the dual representation of U(l), with elements 

{ntn\W^^ = e'"^(nt,| V G [0, 27r). (3.49) 



Wj, as 



Example 1: Consider a two-dimensional space V that decomposes as V = Vo©Vi, where 
the irreps n = and n = 1 are non-degenerate (i.e. do — di — 1). Then the orthogonal 
vectors {\n — 0,to — 1) , \n — 1, ti — 1)} form a basis of V. In column vector notation, 







= |n = 0,to = 1), 







= |n = 1,^1 = 1), 



(3.50) 



the particle number operator n and transformation Wu, read 



n 




1 



1 

e-'"^ 



(3.51) 



Example 2: Consider a four-dimensional space V that decomposes as V = Vq® Vi ® V2, 
where do = d2 = I and di = 2, so that now irrep n = 1 is twofold degenerate. Let 
{|n = 1, ti = 1), |n = 1, ti = 2)} form a basis of Vi. In column vector notation, 
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= \n 


= 0,to = l), 


1 





\n = = 1) 
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= \n 













|n = 2,^2 = 1) 




w 








W 






the particle number operator 


fi and transformation 


read 




/o 


0^ 




fi 





\ 


n 




1 





1 


, W = 











e-'^ 







2j 









e-'2'^ J 



(3.53) 



(3.54) 
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3.2.3.2 Symmetric States and Operators 

In this work we are interested in states and operators that have a simple transformation 
rule under the action of U(l). A pure state G V is symmetric if it transforms as 

W^\^) = e"'"'^!^) V <^ G [0, 2n). (3.55) 

The case n = corresponds to an invariant state, = I^P), which transforms trivially 

under U(l), whereas for n ^ the state is covariant, with being multiplied by a non- 
trivial phase e"^"*^. Notice that a symmetric state |\E') is an eigenstate of n: that is, it 
has a well-defined particle number n. can thus be expanded in terms of a basis of the 
relevant subspace V„, 

n\^) = n\^), 1^) = J2i^n)tjnt„), (3.56) 

tn=l 

where we have introduced a charge label „ on the state coefficients of so that we can 
explicitly associate each coefficient {'^n)tn with its corresponding basis vector |nt„). 

A linear operator T : V — )■ V is invariant if it commutes with the generator h , 

[f,ri] = 0, (3.57) 

or equivalently if it commutes with the action of the group, 

W^fW^^ = f V^G[0,27r). (3.58) 

It follows that T decomposes as (Schur's lemma) 

T = 0f„ (3.59) 

n 

where Tn is a (i„ x dn matrix that acts on the subspace V„ in Eq. (3.44). 

Notice that the operator T in Eq. (3.59) transforms vectors with a well-defined particle 
number n into vectors with the same particle number. That is, U(l)-invariant operators 
conserve particle number. 
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Example 1 revisited: In Example 1 above, symmetric vectors must be proportional to 
either |n = 0,to = l) oi \n = l,ti = 1). An invariant operator T = Tq(BTi is of the form 



T 




, ao, a\ G C. 



(3.60) 



Example 2 revisited: In Example 2 above, a symmetric vector must be of the form 



1^) 








1^) 



/o\ 



or 



/o\ 




V«2/ 



where «i, ^2 ^ C. An invariant operator T = To © Ti © T2 is of the form 



T 











o\ 







/3i 








7i 


5i 













OL2) 



(3.61) 



(3.62) 



where Ti corresponds to the 2x2 central block and ao; cti; A, 7i, ^i, «2 ^ C. 

The above examples illustrate that the symmetry imposes constraints on vectors and oper- 
ators. By using an eigenbasis {|^tn)} of the particle number operator n, these constraints 
imply the presence of the zeros in Eqs. (3.60)-(3.62). Thus, a reduced number of complex 
coefficients is required in order to describe U(l)-symmetric vectors and operators. As we 
will discuss in Sec. 3.2.4, performing manipulations on symmetric tensors can also result 
in a significant reduction in computational costs. 



3.2.3.3 Tensor Product of Two Representations 

Let V(^) and V^^^ be two spaces that carry representations of U(l), as generated by 
particle number operators h^^'' and •h^^\ and let 



(3.63) 



riB 
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be their decompositions as a direct sum of (possibly degenerate) irreps. Let us also 
consider the action of U(l) on the tensor product V^"^^-* = V'-"^^ ® V*-^^ as generated by 
the total particle number operator 

^(AB) (3.64) 
that is, implemented by unitary transformations 

The space V(^^) also decomposes as the direct sum of (possibly degenerate) irreps, 

"AS 

Here the subspace 'Vi^\ with total particle number tiab, corresponds to the direct sum 
of all products of subspaces and YifJ such that + = uab, 

For each subspace yI^^ in Eq. (3.66) we introduce a coupled basis {InABtn^g)}, 

n^'^^^luABtnAB) = nAB\nABtnAB) y {3.68) 

where each vector \nABtnAB) corresponds to the tensor product \nAtnAi ''^Btng) = I^^A^n^)® 
I'^B^ns) of a unique pair of vectors \nAtnA) In^tn^)' with nA + ub = tiab- Let table 
■y-iusc^ with components 

'^n7t„A,nBtr.B^nABtr.AB = AbKabW AtnA^ ^ bKb) , (3.69) 

encode this one-to-one correspondence. Notice that each component of T'"'° is either a 
or a 1. Then 

\nABtnAB) = X! ^"Ttn^.nstns^n^stn^s I^^WS^B^nis)- (3-70) 

For later reference (Sec. 3.4), we notice that T^"'"' can be decomposed into two pieces. 
The first piece expresses a basis {jn^tn^; nBtng)} ofY^"^^^ in terms of the basis {|n^t„^)} 
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of V(^) and the basis {Inst^^)} of V^^^ This assimment occurs as in the absence of the 
symmetry, where one creates a composed index d = b x c hj running, for example, fast 
over index c and slowly over index b as in Eq. (3.26). Note that this procedure does not 
always lead to the set {IriA^nA! '^B^ns)} being ordered such that states corresponding to 
the same total particle number hab = ua + ub are adjacent to each other within the 
set. This ordering is achieved by the second piece: a permutation of basis elements that 
reorganizes them according to their total particle number hab, so that they are identified 
in an one-to-one correspondence with the coupled states {\nABtnAB)}- 

Finally, the product basis can be expressed in terms of the coupled basis 



nAtnAl^Btng) — ^ nABtnAB^"'AtnA,-n-Btng I'^AB^nAs)) (3-71) 
^ABiuAB 



with 

'V^ split ^ fuse / O 

UABtriAB'^^AtnA^'^Btng ~ riAtn a '"-Btng ^riABtuAB ' V'^-'^/ 



Example 3: Consider the case where both Y^^^ and V^^-* correspond to the space of 
Example 1, that is V(^) ^ v[,^^ ©V^^^ and V^^) ^ v{,^^ ®y[^\ where v[,^\ y['^\ \ and 
v[^^ all have dimension 1. Then V^'^'^-' corresponds to the space in Example 2, namely 



where 



= (v(^^©vs-^))®(v;,^)©vf)) 

= Vj,^''^ © VS^""^ © V^^""^ , (3.73) 



Vi^^^ ^ V[,^^ © Vj,^^ (3.74) 

y{AB) ^^(A) ^^(B)^ (3.76) 
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The coupled basis {I^AB^nAs)} reads 



nAB 


= 0,to 


= 1) = 


\nA 


= 0,to 


= 1) ® 




= 0,to 


= 1) 


nAB 


= 


= 1) = 


\nA 


= 0,to 


= 1) 8) 




= 


= 1) 


nAB 


= l,ti 


= 2) = 


\nA 


= 


= 1) ® 


\nB 


= 0,to 


= 1) 


nAB 


= 2,t2 


= 1) = 


\nA 


= 


= 1) ® 


\nB 


= 


= 1), 



where we emphasize that the degeneracy index takes two possible values for uab = 1, 
i.e. ti G {1,2}, since there are two states |nyitn^) ® Instng) with ua + nB = 1. The 
components ^ n^tA,nBtB-*nABtAB ^'^^ tensor T'""" that encodes this change of basis are 
all zero except for 

fuse fuse fuse 'Y^ fuse 1 

-■-Ol.Ol-^Ol ~ 01,11-^11 ~ ^11,01^12 ~ ^11,11^21 ~ ^• 



3.2.3.4 Lattice Models With U(l) Symmetry 

The action of U(l) on the threefold tensor product 

as generated by the total particle number operator 
induces a decomposition 

nABC 

in terms of irreps V^^^^^ which we can now relate to , yH^J and yIt^J . For example, 
we can consider first the product vkl^^ = eg) Y^^J and then the product Y^h]^^^ = 

(AH) (C) 

YriAB ® Ync 5 using a different table T at each step to relate the coupled basis to the 
product basis as discussed in the previous Section. Similarly we could consider the action 
of U(l) on four tensor products, and so on. 

In particular we will be interested in a lattice C made of L sites with vector space Y'^^, 
where for simplicity we will assume that each site s G £ is described by the same finite- 
dimensional vector space V (see Sec. 3.2.2.5). Given a particle number operator n defined 
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on each site, we can consider the action of U(l) generated by the total particle number 
operator 

L 

N = ^&\ (3.81) 
which corresponds to unitary transformations 

Wf^ = e-'^'^ = (e-'""^)®^ = (W^) . (3.82) 
The tensor product space V®^ decomposes as 

V®i^0V^ (3.83) 

N 

and we denote by {\NtN)} the particle number basis in V®^. 

We say that a lattice model is U(l)-symmetric if its Hamiltonian H : Y ^ Y commutes 
with the action of the group. That is, 

[H,N] = 0, (3.84) 

or equivalently 

(W^^ H(W^^^ =H V<^e[0,27r). (3.85) 

One example of a U(l)-symmetric model is the hard core Bose-Hubbard model, with 
Hamiltonian 

L L 

-^HCBH = ^ (^altts+i + Osol+i + ^nsfis+i^ - fx'^Tis, (3.86) 

s=l s=l 

where we consider periodic boundary conditions (by identifying sites L+1 and 1), and a J 
and cis are hard-core bosonic creation and annihilation operators, respectively. In terms 
of the basis introduced in Example 1 these operators are defined as 

/o A /o o' 

a = , n = a'a — \ 

\o oy \ 1 

To see that ^hcbh commutes with the action of the group, we first observe that for two 
sites 

a{a2 + alai , fii + = 0, (3.87) 
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from which it readily follows that 



-f^HCBH, N 



0. 



Notice that the chemical potential term — /i fis = —fiN also commutes with the rest 
of the Hamiltonian. The ground state of Hrcbh in a particular subspace or 

particle number sector can be turned into the absolute ground state by tuning the chemical 
potential fi. This fact can be used to find the ground state {"^ff) of any particle number 
sector through an algorithm which can only minimize the expectation value of i^ncBH- 
However, we will later see that the use of symmetric tensors in the context of tensor 
network states will allow us to directly minimize the expectation value of Hucbr in a 
given particle number sector by restricting the search to states 

1^^) = 5^(^^)i^|iVt^) (3.88) 

tjV = l 

with the desired particle number A^. 

Finally, by making the identifications 

n = , a = 

2 ' 2 

where ax, cry, a^ are the Pauli matrices 



























i. Oi 






ax= ^j, cTy= ^ ^j, a, = I I , (3.89) 

one can map ^hcbh to the spin-| XXZ quantum spin chain 

Hxxz ^ Yl (^^^^^''"^ + ^y^^t'^ + Aa^a^^)) , (3.90) 

where we have ignored terms proportional to N and set A = 7/4. In particular, for A = 
we obtain the quantum XX spin chain 

Hxx^t.i^^^cr^^^'^+alf^al;^^), (3.91) 

s=l 

and for A = 1, the quantum Heisenberg spin chain 

Hxxx - E {ai^'^'^ + + aiM'^'^) . (3.92) 

s=l 

In Sec. 3.2.5, the quantum spin models (3.91) and (3.92) will be used to benchmark the 
performance increase resulting from the use of symmetries in tensor networks algorithms. 
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3.2.4 Tensor Networks With U(l) Symmetry 

In this Section we will consider U(l)-symmetric tensors and tensor networks. I will explain 
how to decompose U(l)-symmetric tensors in a compact, canonical form that exploits 
their symmetry, and then discuss how to adapt the set V of primitives for tensor network 
manipulations in order to work in this form. We will also analyse how working in the 
canonical form affects computational costs. 

3.2.4.1 U(l)-Symmetric Tensors 

Let T be a rank-A; tensor with components Ti-^i2...ik- in Sec. 3.2.2.6, we regard tensor T 
as a linear map between the vector spaces Vl^^l and V[°"*] (3.39). This implies that each 
index is either an incoming or outgoing index. On each space V'''], associated with index 
ii, we introduce a particle number operator n*^'-* that generates a unitary representation 
of U(l) given by matrices W^^^ = e"'"'''*^, ip e [0,27r). In the following, we use W^^^ * to 
denote the complex conjug ate ofWi^l 

Let us consider the action of U(l) on the space 

V[n] ^ yfe] ... yl^fe] (3.93) 



given by 



where 



, (3.94) 



W^^ * if e / 



X? = { Zn : (3.95) 



" ' if GO 



That is, X^'^ acts differently depending on whether index ii of tensor T is an incoming or 
outgoing index. We then say that tensor T, with components Ti^^i^.^ik^ is U(l)-invariant 
if it is invariant under the transformation of Eq. (3.94), 

^ (■*"),,. (■*?'),. ■ ■ ■ (^f ),,, - 

11,12, ---Ak 

for all ip G [0, 2tt). This is depicted in Fig. 3.13. 
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Figure 3.13 - (i) Constraint fulfilled by a U(l)-invariant vector. The only allowed particle 
number on the single index is n = 0. (ii) Constraint fulfilled by a U(l)-invariant matrix. It fol- 
lows from Schur's lemma that the matrix is block-diagonal in particle number, (iii) Constraint 
fulfilled by a rank-three tensor with one incoming index and two outgoing indices. 



Example 4: A U(l)-mvariant vector |\E') — that is, a vector with ri|\E') = and com- 
ponents (\l/„=o)to in the subspace Yn=o which corresponds to vanishing particle number 
n = [cf. Eq. (3.56)]— fulfills 



{^n=o)to' = T i^^) (^"=o)*o V ^ G [0, 27r), (3.97) 

to 

in accordance with Eq. (3.55), as shown in Fig. 3.13(i). 
Example 5: A U(l)-invariant matrix T (3.59) fulfills 

= j:{<)J%-)j. (3.98) 



= Tat(wj) y^e[0,2n), (3.99) 

\ / a'a \ / bb' 

a,b 

in accordance with Eq. (3.58) [see Fig. 3.13(ii)]. 

Example 6: Tensor T in Eq. (3.40), with components Tabc where a and b are outgoing 
indices and c is an incoming index, is U(l)-invariant iff 



Ta'b'c' — 

a,b,c 



5:(^«'^')..(*^1..(*^"-)./- p-^") 

a,b,c 

Y.(^^,^(^':?%,t.(w>p% (3,101) 



for all if G [0,27r) [see Fig. 3.13(iii)]. 
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1^) 



(0 



o 



9^^ 



(ii) 



Figure 3.14 - (i) U(l)-covariant vector l^*), with some non-vanishing particle number n / 0. 
Under the action of U(l) on its index, the covariant vector l^*) acquires a phase e"'""^ (3.103). 
(ii) The U(l)-covariant vector |^), with components (^'n)t„, can be represented by a U(l)- 
invariant matrix T with components Tj^j = (^n)ii, where i is a trivial index (|i| = 1) with 
charge n and is decorated by the opposite arrow to ii. 



Further, we say that a tensor Q with components Qi^i2...if, is U(l)-covariant if under the 
transformation of Eq. (3.94) it acquires a non-trivial phase e"^""^, 

Jl,l2,...,«fc 

for all e [0,27r). 

Example 7: A U(l)-covariant vector |\1') — that is, one which satisfies ■h\'$) = n\'$) for 
some n 7^ 0, and has nonzero components {'^n)t„ only in the relevant subspace V„ [cf. 
Eq. (3.56)]— fulfills 

(^n)i„ = e-'"^(^Oi^ V^G[0,27r), (3.103) 
in accordance with Eq. (3.55). (See also Fig. 3.14.) 

Notice that we can describe the rank-A; covariant tensor Q above by a Ta.nk-{k+l) invariant 
tensor T with components 

Tiii2...ii^i — Qjii2...ifc 1^1 (3.104) 

This is built from Q by adding an extra incoming index i, where index i has fixed particle 
number n and no degeneracy (i.e., i is associated to a trivial space = C). We 
refer to both invariant and covariant tensors as symmetric tensors. By using the above 
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construction, in this work we will represent all U(l)-symmetric tensors by means of U(l)- 
invariant tensors. In particular, we represent the non-trivial components {'^n)tn of 
covariant vector in Eqs. (3.55)-(3.56) as an invariant matrix T of size x 1 with 
components Tt„i = (^n)i„- Consequently, from now on we will mostly consider only 
invariant tensors. 

3.2.4.2 Canonical Form For U(l)-Invariant Tensors 

Let us now write a tensor T in a particle number basis on each factor space in Eq. (3.93). 
That is, each index ii, i2, . . ., ik is decomposed into a particle number index n and a 
degeneracy index t„, ii = (ni,t„J, ^2 = (^2,i^n2), • • •, 4 = (^fc,^nj, and 

'^hh—ik — \ Tnin2...nk j ■ (3.105) 

Here, for each set of particle numbers ni,n2, . . . we regard r„j„2 as a tensor with 
components (^'nin2...nfc) • L^t Ni^ and No^t denote the sum of particle numbers 

corresponding to incoming and outgoing indices. 



ni&I ni&O 

The condition for a non- vanishing tensor of the form T^^nj-.-rifc be invariant under U(l), 
Eq. (3.94), is simply that the sum of incoming particle numbers equals the sum of outgoing 
particle numbers. Therefore, a U(l)-invariant tensor T satisfies 

^= ® Tnin2...njNi,,,N,u,- (3.107) 
ni,n2,...,rn. 

[We use the direct sum symbol to denote that the different tensors T„^„2...nj. are sup- 
ported on orthonormal subspaces of the tensor product space of Eq. (3.93).] In compo- 
nents, the above expression reads 



Ti^i2...ik = \Tnxn2...riA ^Nin,Nouf (3.108) 



Here, S^i^^Nout implements particle number conservation: if A^in 7^ Nout, then all com- 
ponents of T„j„2...nfc must vanish. This generalizes the block structure of U(l)-invariant 
matrices in Eq. (3.59) (where T„„ is denoted T„) to tensors of arbitrary rank k. The 
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canonical decomposition in Eq. (3.107) is important, in that it allows us to identify the 
degrees of freedom of tensor T that are not determined by the symmetry. Expressing 
tensor T in terms of the tensors T„i„2...rife with A^in = A^^out ensures that we store T in the 
most compact way possible. 

Notice that the canonical form of Eq. (3.107) is a particular case of the canonical form 
presented in Eq. (3.18) of Sec. 3.1.2 for more general (possibly non- Abelian) symmetry 
groups. There, a symmetric tensor was decomposed into degeneracy tensors [analogous to 
tensors T„^„2,,,„j, in Eq. (3.107)] and structural tensors [generalizing the term S^.^^n^^^ in 
Eq. (3.107)] which can in general be expanded as a trivalent network of Clebsch-Gordan 
(or coupling) coefficients of the symmetry group. In the case of non-Abelian groups, 
where some irreps have dimension larger than 1, the structural tensors are highly non- 
trivial. However, for the group U(l) discussed in this Section (as for any other Abelian 
group) all irreps are one dimensional and the structural tensors are always reduced to a 
simple expression such as SNin,Nout Eq. (3.107). (Nevertheless, in Sec. 3.4 we will resort 
to a more elaborate decomposition of the structural tensors in order to better exploit 
the presence of symmetry in those tensor network algorithms based on iterating a fixed 
sequence of manipulations.) 

3.2.4.3 U(l)-Symmetric Tensor Networks 

In Sec. 3.2.2.6 we saw that a tensor network J\f where each line has a direction (represented 
with an arrow) can be interpreted as a collection of linear maps composed into a single 
linear map T of which is a tensor network decomposition. By introducing a particle 
number operator on the vector space associated to each line of A/", we can define a unitary 
representation of U(l) on each index of each tensor in JV. Then we say that A/" is a U(l)- 
invariant tensor network if all its tensors are U(l)-invariant. Notice that, by construction, 
if A^ is a U(l)-invariant tensor network, then the resulting linear map T is also U(l)- 
invariant. This is illustrated in Fig. 3.15. 

More generally, we can consider a U(l)-symmetric tensor network, made of tensors that 
are U(l)-symmetric (that is, either invariant or covariant). Recall, however, that any 
covariant tensor can be represented as an invariant tensor by adding an extra index 
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Figure 3.15 - A tensor network J\f made of U(l)-invariant tensors represents a U(l)-invariant 
tensor T . This is seen by means of two equalities. The first equality is obtained by inserting 
resolutions of the identity I = W^W^, on each index connecting two tensors in N. The second 
equality follows from the fact that each tensor in M is U(l)-invariant. 

(3.104). Therefore without loss of generality we can restrict our attention to invariant 
tensor networks. 

3.2.4.4 Tensor Network States and Algorithms With U(l) Symmetry 

As discussed in Sec. 3.2.2.5, a tensor network M can be used to describe certain pure 
states G V®^ of a lattice C If A/" is a U(l)-symmetric tensor network then it will 
describe a pure state that has a well-defined total particle number A^. That is, a 
U(l)-symmetric pure state 

iV|^) = A^l^), e-^^'^l^) = e-^^^l^). (3.109) 

In this way we can obtain a more refined version of popular tensor network states such as 
MPS, TTN, MERA, PEPS, etc. As a variational Ansatz, a symmetric tensor network state 
is more constrained than a regular tensor network state, and consequently it can represent 
less states G N®^ . However, it also depends on fewer parameters. This implies a more 
economical description, as well as the possibility of reducing computational costs during 
its manipulation. 

The rest of this Section is devoted to explaining how one can achieve a reduction in 
computational costs. This is based on storing and manipulating U(l)-invariant tensors 
expressed in the canonical form of Eqs. (3.107)-(3.108). We next explain how to adapt 
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the set V of four primitive operations for the tensor network manipulations discussed in 
Sect 3.2.2.4, namely permutation and reshaping of indices, matrix multiplication, and 
factorization. 



3.2.4.5 Permutation of Indices 

Given a U(l)-invariant tensor T expressed in the canonical form of Eqs. (3.107)-(3.108), 
permuting two of its indices is straightforward. It is achieved by swapping the position 
of the two particle numbers of T„^„2...„j, involved, and also the corresponding degeneracy 
indices. For instance, if the rank-3 tensor T of Eq. (3.40) is U(l)-invariant and has 
components 

Tabc {Tn^nBTic I '^n^+7is,nc (3.110) 

when expressed in the particles number basis a = (nyi,t„^), b = {nB,tng), c = {nc,tna), 
then tensor T' of Eq. (3.25), obtained from T by permuting the last two indices, has 
components 

{T )acb — y^nAricnB) ^nA+nB,nc (3.111) 

where 

{^nAncnB) = ( ^"Ansnc ) • (3.112) 



A''"C "B '^UA'-'iB'-^C 



Notice that since we only need to permute the components of those T„^„^„p such that nA+ 
ns = nc, implementing the permutation of indices requires less computation time than 
a regular index permutation. This is shown in Fig. 3.16, corresponding to a permutation 
of indices using MATLAB. 



3.2.4.6 Reshaping of Indices 

The indices of a U(l)-invariant tensor can be reshaped (fused or split) in a similar manner 
to those of a regular tensor. However, maintaining the convenient canonical form of 
Eqs. (3.107)-(3.108) requires additional steps. Two adjacent indices can be fused together 
using the table T*^""" of Eq. (3.69), which is a sparse tensor made of ones and zeros. 
Similarly an index can be split into two adjacent indices by using its inverse, the sparse 
tensor T=p"' of Eq. (3.72). 
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Figure 3.16 - Computation times (in seconds) required to permute and fuse two indices of 
a rank-four tensor T, as a function of the size of the indices. All four indices of T have 
the same size, 5d, and therefore the tensor contains |T[ = 5^(i^ coefficients. The figures 
compare the time required to perform these operations using a regular tensor and a U(l)- 
invariant tensor, where in the second case each index contains five different values of the 
particle number n (each with degeneracy d) and the canonical form of Eqs. (3.107)-(3.108) 
is used. The upper figure shows the time required to permute two indices: For large d, 
exploiting the symmetry of a U(l)-invariant tensor by using the canonical form results in 
shorter computation times. The lower figure shows the time required to fuse two adjacent 
indices. In this case, maintaining the canonical form requires more computation time. Notice 
that in both figures the asymptotic cost scales as 0{d'^), or the size of T, since this is the 
number of coefficients which need to be rearranged. We note that the fixed-cost overheads 
associated with symmetric manipulations could potentially vary substantially with choice of 
programming language, compiler, and machine architecture. The results given here show the 
performance of a matlab implementation of U(l) symmetry. 
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Example 8: Let us consider again the rank-3 tensor T of Eq. (3.40) with components 
given by Eq. (3.110), where a and h are outgoing indices and c is an incoming index. We 
can fuse outgoing index h and incoming index c into an (e.g. incoming) index rf, obtaining 
a new tensor T' with components 

(T')ad=fe,„J 5n„n,, (3.113) 



where n^i = —UB + nc- (The sign in front of comes from the fact that d is an incoming 
index and h is an outgoing index.) The components of T' are in one-to-one correspondence 
with those of T and follow from the transformation 

"^nAno) , , ^ X^nAUBnc) ^ TiBt^-g^nct^ -^notuj^-- (3.114) 



, , - , ----// / / ■^B'^-ngyn-C'^nQ-^njjtnj^ ■ 

where only the case ua = rir, needs to be considered. To complete the example, let us 
assume that the index a is described by the vector space V(^) ^ Vo © Vi © V2 with 
degeneracies do = 1, di = 2, and d2 = 1] index b is described by a vector space V^-^) = 
V_i © Vo without degeneracies, i.e. d^i = do = 1; and index c is described by a vector 
space V(^) ^ Vo © Vi also without degeneracies, d_i = do = l. Then V^^) = V^^) © V^^) 
(and in this example, also V*-^^ = V*-"^-') and Eq. (3.114) amounts to 

f'] = (fr 



111 



^n) = (^i(-i)o, 

/ 22 V / 211 



(^22)^^ - (^2(-l)l 



111 



where we notice that tensor T' is a matrix as in Eq. (3.62). Similarly, we can split incoming 
index d of tensor T' back into outgoing index b and incoming index c of tensor T according 
to 

TnAnBUc) = ^ (-^nAnu), , noinj^-^UBtng ,nct„fj (3.115) 



•I'D if^njj 



which, again, is non-trivial only for —ub + nc = ud and ua + nB = nc- 
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This example illustrates that fusing and splitting indices while maintaining the canonical 
form of Eqs. (3.107)-(3.108) requires more work than reshaping regular indices. Indeed, 
after taking indices b and c into d = bxchj listing all pairs of values bxc, we still need to 
reorganize the resulting basis elements according to their particle number riij. Although 
this can be done by following the simple table given by T'^"^'', it may add significantly to 
the overall computational cost associated with reshaping a tensor. For instance, Fig. 3.16 
shows that fusing indices of invariant tensors can be more expensive than fusing indices 
of regular tensors. 

3.2.4.7 Multiplication of Two Matrices 

By permuting and reshaping the indices of a U(l)-invariant tensor, we can convert it into 
a U(l)-invariant matrix T = 0^^, Tnn'Sny, or simply 

T = 0f;„, (3.116) 

n 

where T„ = T„„. In components, matrix T reads 

{f)ab = (fn) , (3.117) 

where a = (n,tn) and b = {n,t'^). In particular, similar to the discussion in Sec. 3.2.2.2 
for regular tensors, the multiplication of two tensors invariant under the action of U(l) 
can be reduced to the multiplication of two U(l)-invariant matrices. 

Let R and S be two U(l)-invariant matrices, with canonical forms 

^ = 0^„, S = ^Sn. (3.118) 

n n 

Their product T = R ■ S , (3.28), is then another matrix T which is also block diagonal, 

T = 0f„, (3.119) 

n 

such that each block T„ is obtained by multiplying the corresponding blocks Rn and Sn, 

f^ = Rn- ^n- (3.120) 
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Equations (3.116) and (3.120) make evident the potential reduction of computational 
costs that can be achieved by manipulating U(l)-invariant matrices in their canonical 
form. First, a reduction in memory space follows from only having to store the diagonal 
blocks in Eq. (3.116). Second, a reduction in computational time is implied by only having 
to multiply these blocks in Eq. (3.120). This is illustrated in the following example. 

Example 9: Consider a U(l)-invariant matrix T which is a linear map in a space V that 
decomposes into q irreps V„, each of which has the same degeneracy dn = d. That is, 
T is a square matrix of dimensions dq x rfg, with the block- diagonal form of Eq. (3.116). 
Since there are q blocks T„ and each block has size d x d, the U(l)-invariant matrix T 
contains qd'^ coefficients. For comparison, a regular matrix of the same size contains q^d^ 
coefficients, a number greater by a factor of q. 

Let us now consider multiplying two such matrices. We use an algorithm that requires 
0{l^) computational time to multiply two matrices of size I x /. The cost of performing 
q multiplications of d x d blocks in Eq. (3.120) scales as 0{qd^). In contrast the cost of 
multiplying two regular matrices of the same size scales as 0(g^(i'^), requiring q^ times 
more computation time. 

Figure 3.17 shows a comparison of computation times when multiplying two matrices for 
both U(l)-symmetric and regular matrices. 

3.2.4.8 Factorization of a Matrix 

The factorization of a U(l)-invariant matrix T (3.116) can also benefit from the block- 
diagonal structure. Consider, for instance, the singular value decomposition T = USV of 
Eq. (3.35). In this case we can obtain the matrices 

U=^Un ^=0^n V = ^Vn (3.121) 

n n n 

by performing the singular value decomposition of each block T„ independently, 

fn = UnSnVn. (3.122) 

The computational savings are analogous to those described in Example 9 above for 
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Figure 3.17 - Computation times (in seconds) required to multiply two matrices (upper panel) 
and to perform a singular value decomposition (lower panel), as a function of the size of the 
indices. Matrices of size 5d x 5d are considered. The figures compare the time required to 
perform these operations using regular matrices and U(l)-invariant matrices, where for the 
U(l)-invariant matrices each index contains five different values of the particle number n, 
each with degeneracy d, and the canonical form of Eqs. (3.116)-(3.117) is used. That is, each 
matrix decomposes into five blocks of size d x d. For large d, exploiting the block diagonal 
form of U(l)-invariant matrices results in shorter computation time both for multiplication 
and for singular value decomposition. The asymptotic cost scales with d as 0{d'^), while 
the size of the matrices grows as O(d^). We note that the fixed-cost overheads associated 
with symmetric manipulations could potentially vary substantially with choice of programming 
language, compiler, and machine architecture. The results given here show the performance 
of a MATLAB implementation of U(l) symmetry. 
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the multiplication of matrices. Figure 3.17 also shows a comparison of computation times 
required to perform a singular value decomposition on U(l)-invariant and regular matrices 

using MATLAB. 



3.2.4.9 Discussion 

In this Section we have seen that U(l)-invariant tensors can be written in the canonical 
form of Eqs. (3.107)-(3.108), and that this canonical form is of interest because it offers a 
compact description in terms of only those coefficients which are not constrained by the 
symmetry. We have also seen that maintaining the canonical form during tensor manip- 
ulations adds some computational overhead when reshaping (fusing or splitting) indices, 
but reduces computation time when permuting indices (for sufficiently large tensors) and 
when multiplying or factorizing matrices (for sufficiently large matrix sizes). 

The cost of reshaping and permuting indices is proportional to the size |T| of the tensors, 
whereas the cost of multiplying and factorizing matrices is a larger power of the matrix 
size, for example |Tp/^. The use of the canonical form when manipulating large tensors 
therefore frequently results in an overall reduction in computation time, making it a very 
attractive option in the context of tensor network algorithms. This is exemplified in the 
next Section, where we apply the MERA to study the ground state of quantum spin 
models with a U(l) symmetry. 

On the other hand, however, the cost of maintaining invariant tensors in the canonical 
form becomes more relevant when dealing with smaller tensors. In the next Section we 
will also see that in some situations, this additional cost may significantly reduce, or 
even offset, the benefits of using the canonical form. In this event, and in the specific 
context of algorithms where the same tensor manipulations are iterated many times, it 
is possible to significantly decrease the additional cost by precomputing the parts of the 
tensor manipulations that are repeated on each iteration. Precomputation schemes are 
described in more detail in Sec. 3.4. Their performance is illustrated in the next Section. 
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3.2.5 Tensor Network Algorithms With U(l) Symmetry: 
A Practical Example 

In previous Sections we have described a strategy to incorporate a U(l) symmetry into 
tensors, tensor networks, and their manipulations. To further illustrate how the strategy 
works in practice, in this Section we demonstrate its use in the context of the multi-scale 
entanglement renormalization Ansatz, or MERA, and present numerical results from our 
reference implementation of the U(l) symmetry in matlab. 

3.2.5.1 Multi-Scale Entanglement Renormalization Ansatz 

Figure 3.18 shows a MERA that represent states G Y'^^ of a lattice C made of 
L = 18 sites (see Sec. 3.2.2.5). Recall that the MERA is made of layers of isometric 
tensors, known as disentanglers u and isometrics w, that implement a coarse-graining 
transformation. In this particular scheme, isometrics map three sites into one and the 
coarse-graining transformation reduces the L = 18 sites of C into two sites using two 
layers of tensors. A collection of states on these two sites is then encoded in a top tensor 
t, whose upper index o, — 1,2,..., Xtop is used to label Xtop states |^a) e V®^. This 
particular arrangement of tensors corresponds to the 3:1 MERA described in Evcnbly 
and Vidal (2009a). 
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In this Section we will consider a MERA analogous to that of Fig. 3.18 but with Q 
layers of disentanglers and isometries, which we will use to describe states on a lattice C 
made of 2 x 3*^ sites. We will use this variational Ansatz to obtain an approximation to 
the ground state and first excited states of two quantum spin chains that have a global 
internal U(l) symmetry, namely the spin-1/2 quantum XX chain of Eq. (3.91) and the 
spin-1/2 antiferromagnetic quantum Heisenberg chain of Eq. (3.92). Each spin-1/2 degree 
of freedom of the chain is described by a vector space spanned by two orthonormal states 
{| i)) I t)}- Here we will represent them by the states {|0), |1)} corresponding to zero and 
one particles, as in Example 1 of Sec. 3.2.3.1. For computational convenience, we will 
consider a lattice C where each site contains two spins, or states, {| ||), | |t)) I ti); I tt)}- 
Therefore each site of C is described by a space V = Vq © Vi © V2, where = d2 = 1 
and di = 2, as in Example 2 of Sec. 3.2.3.1. Thus, a lattice C made of L sites corresponds 
to a chain of 2L spins. In such a system, the total particle number ranges from 
to 2L. (Equivalently, the z component of the total spin Sz ranges from —L to L, with 
Sz = N- L.) 



3.2.5.2 MERA With U(l) Symmetry 



A U(l)-invariant version of the MERA, or U(l) MERA for short, is obtained by simply 
considering U(l)-invariant versions of all of the isometric tensors, namely the disentanglers 
M, isometries w, and the top tensor i. This requires assigning a particle number operator 
to each index of the MERA. Each open index of the first layer of disentanglers corresponds 
to one site of C. The particle number operator on any such index is therefore given by 
the quantum spin model under consideration. We can characterize the particle number 
operator by two vectors, n and d: a list of the different values the particle number takes 
and the degeneracy associated with each such particle number, respectively. In the case 
of the vector space V for each site of C described above, n = {0, 1, 2} and d = {1, 2, 1}. 
For the open index of the tensor t at the very top of the MERA, the assignment of charges 
is also straightforward. For instance, to find an approximation to the ground state and 
first seven excited states of the quantum spin model with particle number N, we choose 
n = {N} and d = {8}. (In particular, a vanishing corresponds to = L.) 
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Level q 


Particle numbers n 


Degeneracies d 


top 


{N = 54} 


{Xtop} 


3 


{25,26,27,28,29} 


{1,3,5,3,1} 


2 


{7,8,9,10,11} 


{1,3,5,3,1} 


1 


{1,2,3,4,5} 


{1,3,5,3,1} 





{0,1,2} 


{1,2,1} 



Table 3.1 - Example of particle number assignment in a U(l) MERA for L = 54 sites (or 
108 spins). The total bond dimension isx = l + 3 + 5 + 3+ l = 13. The value of Xtop is set 
as described in the text. 

For each of the remaining indices of the MERA, the assignment of the pair (n, d) needs 
careful consideration and a final choice may only be possible after numerically testing 
several options and selecting the one which produces the lowest expectation value of the 
energy. Table 3.1 shows the assignment of particle numbers and degeneracies made to 
represent the ground state and several excited states in a system of L = 2 x 3'^ = 54 sites 
(that is, 108 spins) with total particle number A^ = L = 54 (or 5*^ = 0). Notice that at 
level q of the MERA {q = 1, 2, 3), each index effectively corresponds to a block of riq = 3"^ 
sites of C. Therefore having exactly Uq particles in a block of Uq sites corresponds to a 
density of 1 particle per site of C. The assigned particle numbers of Table 3.1, namely 
[uq — 2,nq — l,nq,nq + l,nq + 2] for level q, then correspond to allowing for fluctuations 
of up to two particles with respect to the average density. The sum of corresponding 
de generacies d — {c^n,-2, c^ng-i, c^nq, c^nq+i, '^ng+2} gives the bond dimension which in 
the example is x = 13. 

In order to find an approximation to the ground state of either Hxx or Hxxx in 
Eqs. (3.91)-(3.92), we set Xtop = 1 and optimize the tensors in the MERA so as to 
minimize the expectation value 

(*|^|*), (3.123) 

where G is the pure state represented by the MERA and H is the relevant 
Hamiltonian. In order to find an approximation to the Xtop > 1 eigenstates of H with 
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Figure 3.19 - Error in ground state energy /S.E as a function of x foi' the XX and Heisenberg 
models with 2L = 108 spins and periodic boundary conditions, in the particle number sector 
N = L {or Sz = 0). The error is calculated with respect to the exact solutions, and is seen 
to decay exponentially with x- 



lowest energies, we optimize the tensors in the MERA so as to minimize the expectation 
value 



Xtop 



(3.124) 



a=l 



The optimization is carried out using the MERA algorithm described in Evcnbly and 
Vidal (2009a), which requires contracting tensor networks (by sequentially multiplying 
pairs of tensors) and performing singular value decompositions. In the present example, 
all of these operations will be performed exploiting the U(l) symmetry. 

Figure 3.19 shows the error in the ground state energy as a function of the bond dimension 
X, for assignments of degeneracies similar to those in Table 3.2. The error is seen to decay 
exponentially with increasing x, indicating increasingly accurate approximations to the 
ground state. 
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X 


Degeneracy 


no. of 


no. of 


ratio 




J* 
a 


coefficients 


coefficients 










(regular) 


(symmetric) 






4 


[0,1,2,1,0] 


1552 


426 


3.6 


1 


8 


[0,2,4,2,0] 


17216 


4714 


3.7 


1 


13 


[1,3,5,3,1] 


115501 


21969 


5.3 


1 


17 


[1,4,7,4, 1] 


335717 


68469 


5.0 


1 


21 


[1,5,9,5,1] 


779965 


166901 


4.7 


1 


30 


[2,7,12,7,2] 


3243076 


639794 


5.1 


1 



Table 3.2 - Number of coefficients required to specify the tensors of a MERA for L = 54 
as a function of the bond dimension x< decomposed according to a degeneracy vector d. A 
comparison is made between regular tensors and U(l)-invariant tensors. 

3.2.5.3 Exploiting the Symmetry 

We now discuss some of the advantages of using the U(l) MERA. 
Selection of Particle Number Sector 

An important advantage of the U(l) MERA is that it exactly preserves the U(l) symmetry. 
In other words, the states resulting from a numerical optimization are exact eigenvectors 
of the total particle number operator (3.81). In addition, the total particle number A^ 
can be pre-selected at the onset of optimization by specifying it in the open index of the 
top tensor i. 

Figure 3.20 shows the energy gap between the ground state and two excited states of an 
XX chain with 2L spins (or L sites), for N = L particles {Sz = 0). One is the first excited 
state which also has N = L particles. The other is the ground state in the sector with 
A^ = L + 1 particles. The two energy gaps are seen to decay with the system size as L~^. 
The ability to pre-select a given particle number A^ means that only two optimizations 
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Figure 3.20 - Decay of energy gaps A with system size L in the XX model. The upper line 
corresponds to the energy gap between the ground state and the first excited state in the 
N = L particle number (or Sz = 0) sector. The lower line corresponds to the energy gap 
A^^+i between the ground states of the N = L and N = L + 1 particle number sectors. 



were required: one MERA optimization for N = L with Xtop = 2 in order to obtain an 
approximation to the ground state and first excited state of Hxx in that particle number 
sector; and one MERA optimization for = L + 1 with Xtop = 1 in order to obtain an 
approximation to the ground state of Hxx in the particle number sector N = L + 1. 

Similar results can be obtained with the regular MERA. For instance, one can obtain 
an approximation to the ground state of a given particle number sector by adding a 
chemical potential term — /i n^'^^ to the Hamiltonian and carefully tuning the chemical 
potential term /i until the expectation value of the particle number N is the desired one. 
However, the regular MERA cannot guarantee that the states obtained in this way are 
exact eigenvectors of A^. Instead the resulting states are likely to have particle number 
fluctuations. 

Figure 3.21 shows the low energy spectrum of the Heisenberg model Hxxx for a periodic 
system of L = 54 sites (or 108 spins), including the ground state and several excited states 
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Figure 3.21 - Low energy spectrum of Hxxx with L = 54 sites (=108 spins). Depicted 
states have spins of zero (x, blue loops), one (+, red loops), or two (o, green loop), and total 
number of particles {N) between 52 and 56. Note that the second and third spin-1 triplets 
are twofold degenerate. 



both in the particle sector = 54 (or 5*2 = 0), and in neighboring particle sectors. Recall 
that Hxxx is actually invariant under a global internal SU(2) symmetry, of which particle 
number is a U(l) subgroup. Correspondingly the spectrum is organized according to irreps 
of SU(2), namely singlets (total spin 0), triplets (total spin 1), quintuplets (total spin 2), 
etc. Again, using the U(l) MERA, the five particle number sectors A^ = 52, 53, 54, 55, 
and 56 can be addressed with independent computations. This implies, for instance, that 
in order to find the gap between the first and fourth singlets, we can simply set A^ = 54 
and Xtop = 9 on the open index of the top tensor t, to accommodate the first four spin-0 
states and five spin-1 states in the A^ = 54 sector, as seen in Fig. 3.21. In order to capture 
the fourth singlet using the regular MERA, we would need to consider at least Xtop = 19 
(at a larger computational cost and possibly lower accuracy), since this state has only the 
19*^^ lowest energy overall. 
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Reduction of Computational Costs 

The use of U(l)-invariant tensors in the MERA also results in a reduction of computational 
costs. 

First, U(l)-invariant tensors, when written in the canonical form of Eqs. (3.107)-(3.108), 
are block-diagonal and therefore require less storage space. Table 3.2 compares the number 
of MERA coefficients that need to be stored in the regular and symmetric case, for different 
choices of particle number assignments relevant to the present examples. 

Second, the computation time required to manipulate tensors is also reduced when using 
U(l)-invariant tensors in the canonical form. Figure 3.22 shows the computation time 
required for one iteration of the energy minimization algorithm of Evenbly and Vidal 
(2009a) (during which all tensors in the MERA are updated once), as a function of the 
total bond dimension x- The plot compares the time required using regular tensors and 
U(l)-invariant tensors. For U(l)-invariant tensors, we display the time per iteration for 
three different levels of precomputation, as described in Sec. 3.4. The figure shows that 
for sufficiently large Xi using U(l)-invariant tensors leads to a shorter time per iteration 
of the optimization algorithm. 

In our reference implementation (written in matlab), using the symmetry without pre- 
computation is seen to only reduce the computation time by about a factor of 2 for the 
largest x under consideration. This is because maintaining the canonical form for U(l)- 
invariant tensors still imposes a significant overhead for the values of x considered. In 
contrast, when using precomputation we obtained times shorter by a factor of 10 or more. 

The magnitude of the overhead imposed by maintaining the canonical form will depend 
on factors such as programming language and machine architecture, but in general more 
significant gains can be obtained by making full use of precomputation. This option, 
however, requires a significant amount of additional memory (see Sec. 3.4), and a more 
convenient middle ground can be obtained by using a partial precomputation scheme. 
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Figure 3.22 - Computation time (in seconds) for one iteration of the MERA energy mini- 
mization algorithm, as a function of the bond dimension x- For sufficiently large x. exploiting 
the U(l) symmetry leads to reductions in computation time. The horizontal line on this graph 
shows that this reduction in computation time equates to the ability to evaluate MERAs with 
a higher bond dimension x- For the same cost per iteration incurred when optimising a stan- 
dard MERA in matlab with bond dimension x = 20, one may choose instead to optimise a 
U(l)-symmetric MERA with partial precomputation and x = 24, or with full precomputation 
and X = 28. 
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3.3 Conclusions 

In this Chapter we have provided a detailed explanation of how a global internal Abelian 
symmetry may be incorporated into any tensor network algorithm. We considered tensor 
networks constructed from tensors which were invariant under the action of the internal 
symmetry, and showed how each tensor may be decomposed according to a canonical form 
into degeneracy tensors (which contain all the degrees of freedom that are not affected by 
the symmetry) and structural tensors (which are completely determined by the symme- 
try). We then introduced a set of primitive operations V which may be used to carry out 
tensor network algorithms using Ansatze such as MPS, PEPS, and MERA, and showed 
how each of these operations can be implemented in such a way that the canonical form 
is both preserved and exploited for computational gain. 

We then demonstrated the implementation of this decomposition for tensors with an inter- 
nal U(l) symmetry, and computed multiple benchmarks demonstrating the computational 
costs and speed-ups inherent in this approach. We found that although maintaining the 
canonical form imposed additional costs when combining or splitting tensor indices, for 
simulations of a sufficiently large scale these costs can be offset by the gains made when 
performing permutations, matrix multiplications, and matrix decompositions. 

Finally, we implemented the MERA on a quantum spin chain with U(l) symmetry. We 
showed that exploitation of this symmetry can lead to a decrease in the computational cost 
by a factor of 10 or more. These gains may be used either to reduce overall computation 
time or to permit substantial increases in the MERA bond dimension x, and consequently 
in the accuracy of the results obtained. 

Although in this Chapter we have focused on an example which is a continuous Abelian 
group, the formalism presented here may equally well be applied to a finite Abelian group. 
In particular let us consider a cyclic group Zq, q G Z"^.^*^ As in the case of U(l), the Hilbert 
space decomposes under the action of the group into a direct sum of one dimensional 
irreps which are each characterized by an integer charge a, and consequently most of the 

^^Thc fundamental theorem of Abelian groups states that every finite Abelian group may be expressed 
as a direct sum of cyclic subgroups of prime-power order. 
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analysis presented in this Chapter remains unchanged. In particular, matrices which are 
invariant under the action of the group will be block diagonal in the basis labeled by 
charge according to Eq. (3.59), and symmetric tensors enjoy the canonical decomposition 
stated in Eqs. (3.107)-(3.108). The only objects which need modification are the fusion 
and splitting maps, which need to be altered so that they encode the fusion rules for Zq 
instead for U(l). For a cyclic group Z^, the fusion of two charges a and a' gives rise to a 
charge a" according to a" = (a + a')|g where \q indicates that the addition is performed 
modulo q. For example, Z3 has charges a = 0, 1,2, and the fusion rules for take the 
form a X a' ^ a" where the value of a" is given in the following table: 







a 









1 


2 








1 


2 


a' 1 


1 


2 





2 


2 





1 



More generally, a generic Abelian group will be characterised by a set of charges (oi, 02, 03, . . .). 
When fusing two such sets of charges (ai, 02, 03, . . .) and {a[, a'2, a'^, . . .) , each charge is 
combined with its counterpart a[ according to the fusion rule of the relevant subgroup. 
Once again, this behaviour may be encoded in a single fusion map T'^"^'' and its inverse 
■y split rjj_^g formalism presented in this Chapter is therefore directly applicable to any 
Abelian group. 



3.4 Supplement: Use of Precomputation in Iterative 
Algorithms 

We have seen that the use of the canonical form given in Eqs. (3.107)-(3.108) to represent 
U(l)-invariant tensors can potentially lead to substantial reductions in memory require- 
ments and in calculation time. We also pointed out, however, that there is an additional 
cost in maintaining an invariant tensor in its canonical form, and that this is associated 
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with the reshaping (fusing and/or sphtting) of its indices. In some situations this ad- 
ditional cost may significantly reduce, or even offset, the benefits of using the canonical 
form. 

In this Section we investigate techniques for reducing this additional cost in the context of 
iterative tensor network algorithms. Many of the algorithms discussed in Sec. 3.2.2.5 are 
iterative algorithms, repeating the same sequence of tensor network manipulations many 
times over. Examples include algorithms which compute tensor network approximations 
to the ground state by minimizing the expectation value of the energy or by simulat- 
ing evolution in imaginary time, with each iteration yielding an increasingly accurate 
approximation to the ground state of the system. 

The goal of this Section is to identify calculations which depend only on the symmetry 
group, and are independent of the variational coefficients of such algorithms. Where these 
calculations are repeated in each iteration of the algorithm, we can effectively eliminate 
the associated computational cost by performing them only once, either during or prior 
to the first iteration of the algorithm, and then storing and reusing these precomputed 
results in subsequent iterations. We will illustrate this procedure by considering the 
precomputation of a series of operations applied to a single tensor T. 

To do this, we begin by revisiting the fusion and splitting tables of Sec. 3.2.3.3 and 
introducing a graphical representation of these objects. We then introduce a convenient 
decomposition of a symmetric tensor into a matrix accompanied by multiple fusion and/or 
splitting tensors, and linear maps F that map one such decomposition into another. These 
linear maps are independent of the coefficients of the tensor being reorganized, and con- 
sequently they are precisely the objects which can be precomputed in order to quicken 
an iterative algorithm at the expense of additional memory cost. Finally we describe two 
specific precomputation schemes, differing in what is precomputed and in how the pre- 
computed data are utilized during the execution of the algorithm, in order to illustrate the 
trade-off between the amount of memory needed to store the precomputation data and 
the associated computational speedup which may be obtained. In practice, the nature of 
the specific implementation employed will depend on available computational resources. 
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Figure 3.23 - (i) Graphical representation of the fusion tensor T*^"'" and the splitting tensor 
T''"'''*. (ii) The tensors T*^"""" and T''"'''* are unitary, and thus yield the identity when contracted 
pairwise as shown, (iii) A fusion tensor decomposed into two parts. The first part (indicated 
by a circle with an arrow) performs the tensor product of input irreps, n^t^ x n^t^. The result 
is an index that labels pairs {nAtA,nB'tB)- The second part (indicated by a rectangle) is a 
permutation that associates each pair {uAtA^nBtB) with a unique {nABtnAs)' corresponding 
to a vector in the cou pled basis of V^^^). 



3.4.1 Diagrammatic Notation of Fusing and Splitting Tensors 

In describing how we can precompute repeated manipulations of this tensor T, we will find 
it useful to employ diagrammatic representations of the fusion and splitting tables T*^""" 
and T'p"' introduced in Sec. 3.2.3.3. These tables implement a linear map between a pair 
of indices and their fusion product, and thus can be understood as trivalent tensors having 
two input legs and one output leg (or vice versa) in accordance with Sec. 3.2.2.6. We choose 
to represent them graphically as shown in Fig. 3.23(i), where the arrow within the circle 
always points toward the coupled index. The linear maps T'"'° and T^""'" are unitary, and 
consequently we impose that the tensors of Fig. 3.23(i) must satisfy the identities given 
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in Fig. 3.23(ii), corresponding to unitarity under the action of the conjugation operation 
employed in diagrammatic tensor network notation (vertical reflection of a tensor and the 
complex conjugation of its components, typically denoted "I"). Our notation also reflects 
the property, first noted in Sec. 3.2.3.3, that T'"'"' and T'^"* may be decomposed into two 
pieces [Fig. 3.23(iii)]. For the fusion tensor, we identify the first piece (represented by 
a circle containing an arrow) with the creation of a composed index using the manner 
we would employ in the absence of symmetry (3.26). The second piece, represented by 
the small square, permutes the basis elements of the composed index, reorganizing them 
according to total particle number. The two components of the splitting tensor are then 
uniquely defined by consistency with the process of conjugation for the diagrammatic 
representation of tensors, and with the unitarity condition of Fig. 3.23(ii). 

These requirements have an important consequence. Suppose the first part of T im- 
plements 6 X c — 7- d by iterating rapidly over the values of b and more slowly over the 
values of c, and b lies clockwise of c on the graphical representation of T*"'°. This then 
means that on the graphical representation of T^^"* which implements (i — > 6 x c, index 
b must lie counterclockwise of c. It is therefore vitally important to distinguish between 
the splitting tensor and a rotated depiction of the fusing tensor. To this end we require 
that when using this diagrammatic notation, all tensors (with the exception of the fusion 
and splitting tensors) must be drawn with only downward-going legs, as seen for example 
in Fig. 3.24, though the legs are still free to carry either incoming or outgoing arrows as 
before. 



3.4.2 Tree Decomposition 

We find it convenient to decompose a rank-/c, U(l)-invariant tensor T, having components 
^iii2...«fc) cis a binary tree tensor network T consisting of a matrix M which we will call 
the root node, and of k — 2 splitting tensors T"''"' as branching internal nodes, with the 
leaf indices of tree T corresponding to the indices {ii,i2, ■ ■ ■ ,ik} of tensor T. We refer 
to decomposition T as a tree decomposition of T. Figure 3.24 shows an example of tree 



3.4 Supplement: Use of Precomputation in Iterative Algorithms 



103 




-6 



Figure 3.24 - Binary tree decomposition of a synnnnetric tensor T having components 
Tj^jjisMisie- The tree T is comprised of a matrix M as the root node, four splitting ten- 
sors as internal nodes, and ii,i2, as its leaf indices. No incoming or outgoing arrows are 
indicated on the indices in the figure, as the decomposition is valid for any such assignment 
of directional arrows. 



decomposition for a rank-6 tensor. It is of the form 



rp _ \ ^ jiV y split y split y split y split /o lorN 

-L iii2i3iii5i6 / j ilia ji-5>ii,j3 j2-s>i4,i6 j3->i2,i3 i4-s-«4,«5 ' \O.LZ.O) 

h >i2 ,i3 ,i4 

where {ji, J25 is, J4} are the internal indices of the tree. 

The same tensor T may be decomposed as a tree in many different ways, corresponding 
to different choices of the fusion tree. As an example we show two different but equiv- 
alent decompositions of a rank-4 tensor in Fig. 3.25. Different choices 7i,72, . . . of tree 
decomposition for tensor T will lead to different matrix representations Mi, M^-, ... of the 
same tensor. Finally, Fig. 3.26 shows how to obtain the tree decompositions from Ti^^^i^i^ 
by introducing an appropriate resolution of the identity, constructed from pairs of fusion 
operators T and splitting operators T^^"* in accordance with Fig. 3.23(ii). 

The representation of a tensor T by means of a tree decomposition is particularly useful 
because many tensor network algorithms may be understood as a sequence of operations 
carried out on tensors reduced to matrix form. For example, consider tensor network algo- 
rithms such as MPS, MERA, and PEPS. When tensors are updated in these algorithms, 
the new tensor is typically created as a matrix, to which operations from the primitive 
set V of Sec. 3.2.2.4 are then applied. When they are decomposed or contracted with 
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Figure 3.25 - Two possible tree decompositions of a rank-4 tensor T. Different choices 
7i,72,--- of tree decomposition for tensor T lead to different matrices Mi,M25-- - foi' the 
same tensor. 

other tensors, this may once again take place with the tensor in matrix form. Any such 
matrix form may always be understood as the matrix component of an appropriate tree 
decomposition T of tensor T, where the sequence of operations reshaping tensor T to 
matrix M corresponds to the contents of the shaded area in Fig. 3.26. 

3.4.3 Mapping Between Tree Decompositions 

Suppose now that we have a tensor T in matrix form Mi, which is associated with a 
particular choice of tree decomposition 7i, and we wish to transform it into another 
matrix form corresponding to another tree decomposition 72- As indicated, this 
process may frequently arise during the application of many common tensor network 
algorithms. The new matrix M2 can be obtained from Mi by means of a series of reshaping 
(splitting/fusing) and permuting operations, as indicated in Fig. 3.27, and this series of 
operations may be understood as defining a map T: 



The map T is a linear map which depends only on the tree structure of 7i and 72, and 
is independent of the coefficients of Mi. Moreover, T is unitary, and it follows from 
the construction of fusing and splitting tensors and the behaviour of permutation of 
indices (which serves to relocate the coefficients of a tensor) that T simply reorganizes the 




(3.126) 
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Figure 3.26 - Tree decompositions of tensor T are obtained by contracting the tensor with 
an appropriate resolution of the identity on its indices, selected according to the desired choice 
of the fusion tree T. In each instance, evaluation of the contents of the shaded region yields 
the appropriate matrix M. 



106 



Abelian Symmetries of Spin Systems 




Figure 3.27 - A matrix Mi can be reorganized into another matrix M2 by means of fusion 
tensors, splitting tensors, and the permutation of indices. These operations define a one to 
one linear map T that acts to reorganize the coefficients of Mi. F does not depend on the 
coefficients of Mi, but solely on the sequence of operations performed. 

coefficients of Mi into the coefficients of M2 in a one-to-one fashion. 

Therefore, one way to compute the matrix M2 from matrix Mi is by first computing the 
hnear map F, which is independent of the specific coefficients in tensor T, and by then 
applying it to Mi. 

3.4.4 Precomputation Schemes For Iterative Tensor Network 
Algorithms 

The observation that the map F is independent of the specific coefficients in Mi is par- 
ticularly useful in the context of iterative tensor network algorithms. It implies that, 
although the coefficients in Mi will change from iteration to iteration, the linear map F 
in Eq. (3.126) remains unchanged. It is therefore possible to calculate the map F once, 
during the first iteration of the simulation, and then to store it in memory and re-use it 
during subsequent iterations. We refer to such a strategy as a precomputation scheme. 
Figure 3.28 contrasts the program flow of a generic iterative tensor network algorithm 
with and without precomputation of the transformations F. 



Using such a precomputation scheme, a signiflcant speed-up of simulations can be ob- 
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Figure 3.28 - Flow diagram for the execution of a predetermined number of iterations of 
a generic iterative tensor network algorithm, (i) without any precomputation and (ii) with 
precomputation of the operations F. 



tained, at the price of storing potentially large amounts of precomputed data (as a single 
iteration of the algorithm may require the application of many different transformations 
r). Therefore a trade-off necessarily exists between the amount of speed-up that can be 
obtained and the memory requirement which this entails. In this Section we describe two 
different precomputation schemes. The first one fully precomputes and stores all maps 
r, and is relatively straightforward to implement. This results in the maximal increase 
in simulation speed, but implementation requires a large amount of memory. The sec- 
ond scheme only partially precomputes the maps F, resulting in a moderate speed-up of 
simulations, but with memory requirements which are also similarly more modest. 

3.4.4.1 Maximal Precomputation Scheme 

As noted in Sec. 3.4.3, applying the map F to a matrix Mi simply reorganizes its coeffi- 
cients to produce the matrix M2. Moreover, if the indices of matrices Mi and M2 are fused 
to yield vectors Vi and V2 then the map F may be understood as a permutation matrix, 
and this in turn may be concisely represented as a string of integers F = 71, . . . ,7|Mi| 
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such that entry i of V2 = FVi is given by entry 7^ of vector Vi. Because all of the ele- 
ments from which F is composed are sparse, unitary, and composed entirely of O's and 
I's, the permutation to which F corresponds may be calculated at a total cost of only 
0(|Mi|), where \Mi\ counts only the elements of Mi which are not fixed to be zero by 
the symmetry constraints of Eq. (3.107). In essence, for each element of the vector Vi 
one identifies the corresponding number and degeneracy indices {nf^^ , tf^'^ ) on each leg 
i G {1,2} of matrix Mi. One can then read down the figure, applying each table T'""" 
or T'p"' in turn to identify the corresponding labels {n',t') on the intermediate legs, until 
finally the corresponding labels on the indices of M2 are obtained. There is then a further 
1:1 mapping from each set of labels {n^^,t1^^), (n^^^,t2^^) on M2 to the corresponding 
entry in V2, completing the definition of F as a map from Vi to V2. 

Storing the map F for a transformation such as the one shown in Fig. 3.27 imposes a 
memory cost of 0(|Mi|). The application of this map also incurs a computational cost of 
0(|Mi|), but computational overhead is saved in not having to reconstruct the map F on 
every iteration of the algorithm. 

3.4.4.2 Partial Precomputation Scheme 

The 0(|Mi|) memory cost incurred in the previous scheme can be significant for large 
matrices. However, we may reduce this cost by replacing the single permutation F em- 
ployed in that scheme with multiple smaller operations which may also be precomputed. 
In this approach Mi is retained in matrix form rather than being reshaped into a vector, 
and we precompute permutations to be performed on its rows and columns. 

First, we decompose all the fusion and splitting tensors into two pieces in accordance with 
Fig. 3.23(iii). Next, we recognise that any permutations applied to one or more legs of 
a fusion or splitting tensor may always be written as a single permutation applied to the 
coupled index [Fig. 3.29(i)]. We use this to replace all permutations on the intermediate 
indices of the diagram with equivalent permutations acting only on the indices of Mi and 
the open indices, as shown for a simple example in Fig. 3.29(ii). The residual fusion and 
splitting operations, depicted by just a circle enclosing an arrow, then simply carry out 
fusion and splitting of indices as would be performed in the absence of symmetry (3.26)- 
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Figure 3.29 - (i) Permutations applied to one or more legs of a fusion or splitting tensor can 
be replaced by an appropriate permutation on the coupled index. This process can be used to 
replace all permutations applied on internal indices of a diagram such as Fig. 3.27 with net 
permutations on the indices of A/i and on the open indices of the network, as in shown in (ii). 
The residual fusion and splitting operations, depicted as an arrow in a circle, simply perform 
the basic tensor product operation and its inverse, (3.26)-(3.27), as described in Fig. 3.23(iii) 
and Sec. 3.4.1. 



(3.27). These operations are typically far faster than their symmetric counterparts as 
they do not need to sort the entries of their output indices according to particle number. 

In subsequent iterations, the matrix M2 is obtained from Mi by consecutively 

1. permuting the rows and columns of Mi using the precomputed net permutations 
which act on the legs of Mi; 

2. performing any elementary (non-symmetric) splitting, permuting of indices, and 
fusing operations, as described by the grey-shaded region in Fig. 3.29(ii); 

3. permuting the rows and columns of the resulting matrix, using the precomputed net 
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permutations which act on the open legs of Fig. 3.29 (ii). 

When matrix Mi is defined compactly, as in (3.107), so that elements which are identically 
zero by symmetry are not explicitly stored, a tensor T is constructed from multiple blocks 
identified by U(l) charge labels on their indices [T„^„2...„^ in Eq. (3.107)]. Under these 
conditions the elementary splitting, fusing, and permutation operations of step 2 above are 
applied to each individual block, but some additional computational overhead is incurred 
in determining the necessary rearrangements of these blocks arising out of the actions 
performed. This rearrangement may be computed on the fiy, or may also be precomputed 
as a mapping between the arrangement of blocks in Mi and that in M2. 

The memory required to store the precomputation data in this scheme is dominated 
by the size of the net permutations collected on the matrix indices, and is therefore of 
0{^\ Ml I ) . The overall cost of obtaining M2 from Mi is once again of 0(|Mi|), but 
is in general higher than the previous scheme as this cost now involves two complete 
permutations of the matrix coefficients, as well as a reorganisation of the block structure 
of Ml which may possibly be computed at runtime. Nevertheless, in situations where 
memory constraints are significant, partial precomputation schemes of this sort may be 
preferred. 

3.5 Supplement: Notes on the Implementation of 
Abelian Symmetries 

As noted in Sec. 3.1.2, a symmetric tensor may be decomposed into a spin network 
and a collection of degeneracy tensors. This was seen again in Sees. 3.2.3-3.2.4 for the 
symmetry group U(l), with the spin network in this instance being trivial due to the 
Abelian nature of the group. For an Abelian symmetry, we may consequently understand 
this decomposition [(3.59), (3.107)] as dividing a tensor into a number of blocks, labelled 
by the charges on each leg of the tensor, the majority of which are systematically zero. 
An example of this is seen in Eq. (3.62) for a two-legged tensor, and is reproduced here 
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with the different charge blocks highlighted and labelled: 



T 



bjO 
o 



charge 
1 



v 



ao 

ai I3i 

7i 5i 

J 



(3.127) 



For an Abelian tensor, this block decomposition may be trivially extended over an arbi- 
trary number of legs, with each block being indexed by the associated charges on all legs 
(3.105), and this fact may be exploited in implementation of the symmetry. 

To illustrate how this is done, consider the fusing of two indices under the action of the 
simplest Abelian symmetry group, Z2, and recall that fusing and splitting of indices may 
be thought of as consisting of two stages (see Sees. 3.2.3.3 and 3.4.1). In the first stage, 
the charges on the legs are combined according to a typical tensor product process, as per 
Eq. (3.26), iterating rapidly over one charge and slowly over the other: 

+ X + — ^ + 
- X > +. 



(3.128) 



In the second stage, the states of the tensor product space are re-ordered to collect like 
charges: 



+ X 

— X 

+ X 

— X 



+ 



(3.129) 



Significantly, this entire process may be conducted at the level of blocks, and no mixing 
or re-ordering of the elements within individual blocks is required. If the entire tensor is 
maintained in the form of a sparsely-populated array of blocks in this manner, and blocks 
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Figure 3.30 - A Z2-symmetric tensor with three indices undergoes fusion, first of indices 2 
and 3, then of indices 1 and 2 (as indicated by the upper arrows). This is achieved via a 
re-organisation of its block structure, and re-shaping of individual blocks. Blocks with non- 
trivial content in each diagram are shaded. We track the location of one individual block, 
marked *, and note that if in diagram (i) it has dimension xi x 2:2 x X3, then in (ii) it will 
have dimension xi x (X2X3), and in (iii), dimension (X1X2X3) x 1. The lower arrows indicate 
the reverse process. 



are not concatenated on fusion, then both fusion and sphtting may be performed without 
requiring the addressing of the contents of any individual blocks. This is illustrated 
in Fig. 3.30. [It is assumed either that any reshaping of individual blocks is deferred, in 
keeping with the philosophy of Sec. 3.4, or that an efficient representation of n-dimensional 
tensors employed — such as that used by matlab — for which reshape operations may be 
performed on the individual blocks at trivial cost independent of the size of the block.] A 
similar treatment may be applied to all the primitives of set V (Sec. 3.2.2.4). 

Using this approach, the author was able to implement the Z2-symmetric ternary ID 
MERA with a demonstrable increase in performance over the standard MERA for bond 
dimensions x > 12. Performance for Abelian symmetries with greater numbers of charge 
sectors, such as U(l), is anticipated to be even higher, effectively removing the need for 
the precomputation techniques discussed in Sec. 3.4. However, these techniques remain 
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important for the exploitation of non- Abelian symmetries (Chapter 5 and Singh and Vidal 
in preparation) and anyons (Chapter 4 and Pfeifer et al. 2010). 

One note of caution: When combining indices using a symmetry-preserving ordering such 
as that given in Eq. (3.129), care must be taken to ensure that fusion is implemented in 
an associative manner, i.e. d = {{a x b) x c) results in the same ordering of entries on 
index d as the fusion c? = (a x (6 x c)), if consistent results are to be obtained. 

3.5.1 Fermions 

My initial development of efficient Z2 symmetry algorithms was performed in parallel with 
Phihppe Corboz, who adopted a similar computational strategy. To avoid inappropriate 
duplication of research efforts, it was decided that Philippe would then proceed to incor- 
porate fermionic exchange behaviour, while I would instead study the extension of this 
technique to other symmetry groups, resulting in the present work on Abelian symme- 
tries, non-Abelian symmetries, and anyons. The principle by which fermionic exchange 
statistics may be incorporated into this scheme is, however, easily understood as follows: 

Recall that we have chosen to enumerate the legs of a tensor counterclockwise from the 
9 o'clock position (Sec. 3.2.2.1). When a tensor network diagram involves crossings of pairs 
of legs, these crossings may be understood as a permutation operation which exchanges 
the ordering of the indices on the tensor. However, for fermionic statistics, we have a 
Z2 symmetry and thus for any given index, each block in the tensor will have a charge 
label, + or — , associated with that index. The permutation operation corresponding to 
index exchange is performed independently on each block in the normal fashion for a 
Z2-symmetric tensor, subject to one simple modification: When the charge labels for a 
given block are — on both indices, in addition to the permutation operation, the block is 
multiplied by —1. This is the essence of the "swap gate" formalism introduced by Corboz 
et al. (2010b)— see also Barthel et al. (2009); Corboz and Vidal (2009); Pineda et al. 
(2010); Pizorn and Verstraete (2010). 

I note that even greater efficiency may be attained not only by writing each tensor as a 
grid of blocks, but also by associating with each block a numeric multiplier, initially 1. To 
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d b c f 

Figure 3.31 - For an Abelian group symmetry, the contraction of tensors R and S illustrated 
in Fig. 3.8 is reduced to the matrix multiplication T" = R" x S" , which may be represented 
as shown in diagram (i), above [or diagram (3) of Fig. 3.8]. However, there is an implicit 
index exchange present in this multiplication, shown explicitly in (ii), and for fermions this will 
introduce additional factors of —1. 



multiply a block by —1, it now suffices to multiply the corresponding numeric factor by 
— 1. In this manner the cost of the swap gate is reduced, requiring only one operation to 
change the sign of an entire block. Of course, this minus sign must eventually be applied 
in order to obtain numerical results, but if many swap gates are applied to a tensor before 
these numerical results are computed, then deferring (and combining) their evaluation in 
this manner may result in a significant saving in calculation time. A generalised version of 
this technique is also used to reduce computational cost in the study of anyonic systems, 
as described in Sec. 4.7 

One caution is required when using the methods of Sec. 3.2 for the simulation of fermions. 
Note that when contracting a pair of tensors as per Sec. 3.2.2.2, the counterclockwise 
ordering of the indices to be contracted must be the same on each tensor. There will 
frequently be index permutations required to set up this ordering, which will introduce 
some exchange factors of —1. However, there are also exchange factors associated with the 
contraction itself, as seen in e.g. Fig. 3.31, and these too must be taken into account. For 
an alternative approach where exchange factors are associated only with the re-ordering 
of indices, and there are no such hidden factors associated with the contraction of two 
tensors, see the discussion in Sec. 5.2. 



Chapter 4 



Anyonic Tensor Networks 



Having developed a formalism for the exploitation of internal symmetries in tensor net- 
works and exploited it for Abelian symmetry groups, it was natural to seek to apply this 
formalism to more general classes of physical systems. As a colleague (Sukhwinder Singh) 
was already working on the implementation of SU(2) symmetries for the MPS and MERA, 
with obvious extension to spin systems having any non-Abelian internal symmetry group, 
it was decided that I would instead study whether our approach could be extended to 
permit the study of systems of anyons. 

Anyonic systems do not in general exhibit an internal symmetry group; instead their 
behaviours and statistics are typically described by a more general structure known as a 
Unitary Braided Tensor Category (UBTC). To specify an anyon model using a UBTC, one 
must declare a set of charge labels (including a vacuum charge), fusion rules describing 
how they combine (which must be associative), a set of basis transformations known as 
F moves, and a tensor which describes the exchange statistics of the charges. These 
properties will all be described in detail in Sec. 4.2. For any group it is possible to 
construct an associated UBTC where the F moves are derived from the 6-j symbols, and 
the tensor reflects the universal braid matrix for the irreps of the group. In fact, 
this UBTC provides a complete description of the group at the level of representations, 
up to (but not including) explicit construction of representations of the irreps themselves 
(see Chapter 5). However, although every group may be associated with a UBTC, the 
converse does not hold, and UBTCs may also be used to describe more general structures 
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such as quantum doubles [e.g. D(D3)], and quantum groups [e.g. g-deformed SU(2), or 
SU(2)fc]. It is the UBTCs associated with these more general structures which may be 
used to describe systems with anyonic statistics. 

Because of the close association between groups and UBTCs, it was tempting to ask 
whether our internal symmetry formalism, which could be applied to quantum systems 
with a group-based mathematical structure, could be extended to address systems based 
on any UBTC Of particular interest are the UBTCs associated with quantum groups, as 
these are associated with many interesting anyon models such as SU(2)3, for which the 
integer subalgebra describes a class of anyons known as "Fibonacci anyons", capable of 
supporting universal quantum computation simply through particle exchange. However, 
many of the charges in these models may not be associated with explicit matrix represen- 
tations, as in many cases the fusion rules imply that the dimensions of these irreps must 
be non-integer, or even irrational. Consequently the problem could not be viewed as one 
of decomposing a known system into charge sectors as per Chapter 3, but instead had to 
be formulated directly in the graphical language of UBTCs which constitutes a natural 
description of an anyonic system. 

In this Chapter, I describe a formalism permitting exactly this, whereby any tensor net- 
work Ansatz or algorithm may be constructed for a system of anyons, or indeed for any 
other model describable in terms of a UBTC. This Chapter will only concern itself ex- 
plicitly with systems of anyons on the disc, though a subsequent treatment of anyons on 
surfaces of higher genus is planned. 

(The same formalism may even be applied to ID systems described by a Unitary Tensor 
Category without including a notion of particle exchange, though the author is as yet 
unaware of any physical systems of interest which take such a form.) 



Sections 4.1-4.6 of this Chapter have previously been published as Pfcifcr, Corboz, Buer- 
schaper, Aguado, Troyer, and Vidal, Physical Review B, 82, 115126, 2010, (c) (2010) by 
the American Physical Society. 
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The study of anyons offers one of the most exciting challenges in contemporary physics. 
Anyons are exotic quasiparticles with non-trivial exchange statistics, which makes them 
difficult to simulate. However, they are of great interest as some species offer the prospect 
of a highly fault-tolerant form of universal quantum computation (Kitacv, 2003; Nayak 
et al., 2008), and it has been suggested that the simplest such species may appear in the 
fractional quantum Hall state with filling fraction v = 12/5 (Xia et al., 2004). Despite the 
current strong interest in the development of practical quantum computing, our ability 
to study the collective behaviour of systems of anyons remains limited. 

The study of interacting systems of anyons using numerical techniques was pioneered by 
Feiguin et al. (2007), using exact diagonalisation for ID systems of up to 37 anyons, and 
the Density Matrix Renormalisation Group algorithm (DMRG) (White, 1992) for longer 
chains. Also related is work by Sierra and Nishino (1997), later extended by Tatsuaki 
(2000), which applies a variant of DMRG to spin chain models having SU{2)k symmetry. 
Some of these models are now known to correspond to SU{2)k anyon chains (Trcbst et al., 
2008a), and using this mapping these systems may also be studied using the Bethe Ansatz 
(Alcaraz et al., 1987) and quantum Monte Carlo (Todo and Kato, 2001). 

However, all of these methods have their limitations. Exact diagonalisation has a com- 
putational cost which is exponential in the number of sites, strongly limiting the size of 
the systems which may be studied. DMRG is capable of studying larger system sizes, but 
is typically limited to ID or quasi-lD systems (e.g. ladders). Mapping to a spin chain 
is useful in one dimension but is substantially less practical in two. There are therefore 
good reasons to desire a formalism which will allow the application of other tensor net- 
work algorithms to systems of anyons. Many of these tensor networks, such as Projected 
Entangled Pair States (PEPS) (Gu et al., 2008; Jordan et al., 2008; Nishino and Oku- 
nishi, 1998; Vcrstractc and Cirac, 2004; Xie et al., 2009), and the 2D versions of Tree 
Tensor Networks (TTN) (Tagliacozzo et al., 2009) and of the Multi-scale Entanglement 
Renormalisation Ansatz (MERA) (Cincio et al., 2008; Evenbly and Vidal, 2009b, 2010a) 
have been designed specifically to accurately describe two-dimensional systems. 
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In one dimension, many previously studied systems of interacting anyons display extended 
critical phases (e.g. Feiguin ct al., 2007; Trebst et al., 2008a), which are characterised by 
correlators exhibiting polynomial decay (Di Francesco et al., 1997). Whereas DMRG 
favours accurate representation of short range correlators at the expense of long-range 
accuracy, the ID MERA (Vidal, 2007a, 2008) is ideally suited to this situation as its 
hierarchical structure naturally encodes the renormalisation group flow at the level of 
operators and wavefunctions (Chen et al., 2010a; Vidal, 2007a, 2008, 2010), and hence 
accurately reproduces correlators across a wide range of length scales (Evenbly and Vidal, 
2009a; Giovannetti et al., 2008; Pfeifer et al., 2009; Vidal, 2007a, 2008). The development 
of a general formalism for anyonic tensor networks is therefore also advantageous for the 
study of ID anyonic systems. 

This Chapter describes how any tensor network algorithm may be adapted to systems of 
anyons in one or two dimensions using structures which explicitly implement the quantum 
group symmetry of the anyon model. As a specific example I demonstrate the construction 
of the anyonic ID MERA, which I then apply to an infinite chain of interacting Fibonacci 
anyons at criticality. The approach which I present is completely general, and can be 
applied to any species of anyons and any tensor network Ansatz. 



Consider a lattice Cq of n sites populated by anyons. In contrast to bosonic and fermionic 
systems, for many anyon models the total Hilbert space can not be divided into a 
tensor product of local Hilbert spaces. Instead, a basis is defined by introducing a specific 
fusion tree [e.g. Fig. 4.1(i)]. The fusion tree is always constructed on a linear ordering of 
anyons, and while the ID lattice naturally exhibits such an ordering, for 2D lattices some 
linear ordering must be imposed. Each line is then labelled with a charge index such 
that the labels are consistent with the fusion rules of the anyon model. 



For anyon types where some entries of the multiplicity tensor N^^^ take values greater than 
1, a label Ui is also affixed to the vertex which represents the fusion process to distinguish 
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Figure 4.1 - (i) Example representation of a state \ijj) in a fusion tree basis for a system of 6 
anyons. Labels ai indicate charges associated with edges of the fusion tree graph, and labels 
Ui are degeneracies associated with vertices. The structure of the tree corresponds to a choice 
of basis, and does not affect the physical content of the theory, (ii) Braiding may be used 
to change the ordering of the leaves of a fusion tree basis, or to represent anyon exchange, 
(iii) F-moves convert between the bases associated with different fusion trees. 



between the different copies of charge c. The edges of the graph which are connected to a 
vertex only at their lower end are termed "leaves" of the fusion tree, and we will associate 
these leaves with the charge labels ai . . . a„. Different orderings of the leaves on a fusion 
tree may be interconverted by means of braiding [Fig. 4.1(ii)], and different fusion trees, 
corresponding to different bases of states, may be interconverted by means of F moves 
[Fig. 4.1 (iii)] (Bonderson, 2007; Kitaev, 2006). In some situations it may also be useful to 
associate a further index hi with each of the leaves of the fusion tree. For example, if the 
leaves are equated with the sites of a physical lattice, then this additional index may be 
used to enumerate additional non-anyonic degrees of freedom associated with that lattice. 
For simplicity we will usually leave these extra indices 6i . . . 6„ implicit, as we have done 
in Fig. 4.1, as they do not directly participate in anyonic manipulations such as F moves 
and braiding. 
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Let the total number of charge labels on the fusion tree be given by m, where m > n. For 
Abelian anyons the fusion rules uniquely constrain all for i > n, and provided there 
are no constraints on the total charge, the total Hilbert space reduces to a product of 
local Hilbert spaces V, such that Yco = V®". For no n- Abelian anyons, additional degrees 
of freedom arise because some fusion rules admit multiple outcomes, permitting certain 
tti {i > n) to take on multiple values while remaining consistent with the fusion rules, and 
the resulting Hilbert space does not necessarily admit a tensor product structure. 

We will now associate a parameter i/j with each charge on the fusion tree, which we 
will term the degeneracy. This parameter corresponds to the number of possible fusion 
processes by which charge may be obtained at location i. Where charge ak arises from 
the fusion of charges a, and aj, then i'k,ak will satisfy 



For systems where the only degrees of freedom are anyonic, degeneracies on the physical 
lattice Cq (i.e. z/j a., 1 < i < n) will take values of or 1 depending on whether a charge 
ai is permitted on lattice site i. Higher values of z/j^a. may be used on the physical lattice 
if there is also a need to represent additional non-anyonic degrees of freedom, enumerated 
by indices 6i . . . 6„. 

Up to this point we have parameterised our Hilbert space in terms of explicit labellings of 
the fusion tree. We now adopt a different approach: Consider an edge i of the fusion tree 
which is not a "leaf" . As well as labelling this edge with a charge a, we may introduce a 
second index /ij, running from 1 to z/j_a-. Each pair of values {aj,/ij} may be associated 
with a unique charge labelling for the portion of the fusion tree from edge i out to the 
leaves, with these labellings being compatible with the fusion rules in the presence of a 
charge of on site i (for an illustration of this, see Fig. 4.2). Provided we know the 
structure of the fusion tree above i and have a systematic means of associating labellings 
of that portion of the tree with values of /ij, then in lieu of stating the values of all aj 
for edges j involved in that portion of the tree, we may simply specify the value of the 
degeneracy index /ij. In this way we may specify an entire state in the form 






(4.3) 
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Figure 4.2 - The leaves of this fusion tree carry the charge labels ai to ag. An edge which 
is not a leaf, labelled with charge aio, is indicated by the large grey arrow. The portion of 
the fusion tree extending from edge aio out to the leaves is indicated by the grey ellipse. If a 
degeneracy index /iio is associated with charge aio, then for a given value of aio, index /_iio 
will enumerate all compatible labellings of the highlighted portion of the fusion tree. 

where is the total charge obtained on fusing all the anyons. The index /i^, which is the 
degeneracy index associated with the total charge of the fusion tree, may be understood 
as systematically enumerating all possible labellings of the entire fusion tree including 
charge labels, vertex labels, and any labels associated with additional non-anyonic degrees 
of freedom. For an example, see Fig. 4.3. Note that for a given edge i, the value of the 
degeneracy z/j^a. may vary with the charge and consequently the range of the degeneracy 
index /i^ in Eq. (4.3) is dependent on the value of the charge a^- 

The notation of Eq. (4.3) should be contrasted with that of Fig. 4.1(i). In the latter, 
the number of indices on c depends upon the number of charge labels on the fusion tree, 
whereas in the former, the tensor describing the state is always indexed by just one pair of 
labels — charge and degeneracy — which will prove advantageous in constructing a tensor 
network formalism for systems of anyons. 

We now choose to restrict our attention to systems having the identity charge. We may 
do this without loss of generality as a state on n lattice sites with a total charge may 
always be equivalently represented by a state on n + 1 lattice sites whose total charge is 
the identity, with a charge on lattice site n + 1. This additional charge annihilates the 
total charge am of sites 1 ... n to give the vacuum. The expression for \ip) then becomes 




(4.4) 



122 



Anyonic Tensor Networks 







T T T T 



T 




T T T T 




Figure 4.3 - Example enumeration of states according to am and /x^ for a fusion tree de- 
scribing four Fibonacci anyons. The Fibonacci anyon model has one non-vacuum charge label 
(r) and one non-trivial fusion rule, r x r — )■ 1 -|- r. Because the charges 1 and r are both 
self-dual, no arrows are required on diagrammatic representations of Fibonacci anyon fusion 
trees. 



where /jim' ranges from 1 to the dimension of the Hilbert space of the system of n sites 
with total charge am- Consequently we may represent the state |^) of a system of anyons 
by means of the vector Ci^^, . For simplicity of notation, we will take greek indices from 
the beginning of the alphabet to correspond to pairs of indices {ai,/Xj} consisting of a 
charge index and the associated degeneracy index. The vector Ci^^, will therefore be 
denoted simply c", with the understanding that in this case the charge component a^/ 
of multi- index a takes only the value 1. (Multi-index a is raised as we will shortly 
introduce a diagrammatic formalism in which vector c is represented by an object with 
a single upward-going leg. In this formalism, upward- and downward-going legs may be 
associated with upper and lower multi- indices respectively.) 



4.3 Anyonic Operators 

We will divide our consideration of anyonic operators into two parts. First we shall 
consider operators which map a state on some Hilbert space T-L into another state on the 
same Hilbert space. When apphed to a state represented by c", such an operator leaves 



4.3 Anyonic Operators 



123 



the degeneracies of the charges in multi-index a unchanged. We will therefore call these 
degeneracy-preserving anyonic operators. Then we will consider those operators which 
map a state on some Hilbert space Ti into a state on some other Hilbert space H'. These 
operators may represent processes which modify the environment, for example by adding 
or removing lattice sites, and also play an important part in anyonic tensor networks, 
for instance taking the role of isometrics in the TTN and MERA. As these operators 
can change the degeneracies of charges in a multi-index a, we will call them degeneracy- 
changing anyonic operators. More generally, the degeneracy-preserving anyonic operators 
may be considered a subclass of the degeneracy- changing anyonic operators for which 

n = n'. 

4.3.1 Degeneracy-Preserving Anyonic Operators 

We begin with those operators which map states on some Hilbert space H into other 
states on the same Hilbert space H. Examples of these operators include Hamiltonians, 
reduced density matrices, and unitary transformations such as the disentanglers of the 
MERA. 

First, we introduce splitting trees. The space of splitting trees is dual to the space of 
fusion trees. While the space of fusion trees consists of labelled directed graphs whose 
number of branches increases monotonically when read from bottom to top, the space 
of splitting trees consists of labelled directed graphs whose number of branches increases 
monotonically when read from top to bottom. An inner product is defined by connecting 
the leaves of fusion and splitting trees which have equivalent linear orderings of the leaves 
(braiding first if necessary), then eliminating all loops as per Fig. 4.4(i), with F moves 
performed as required. 

Anyonic operators may always be written as a sum over fusion and splitting trees, such 
as the two-site operator M shown in Fig. 4.4(ii), and for degeneracy-preserving anyonic 
operators it is always possible to choose the splitting tree to be the adjoint of the fusion 
tree. To apply an operator to a state the two corresponding diagrammatic representations 
are connected as shown in Fig. 4.4(iii), and closed loops may be eliminated as shown in 
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(i) 




1 



Figure 4.4 - (i) Loops are eliminated by replacing them with an equivalent numerical fac- 
tor determined by the normalisation convention. The factor given here corresponds to the 
diagrammatic isotopy convention employed in Bonderson (2007). (ii) Definition of a simple 
two-site anyonic operator, (iii) Application of an operator to a state is performed by con- 
necting the diagrams' free legs. By performing F moves and eliminating loops (and in more 
complex examples, also braiding) it is possible to obtain an expression for the resulting state 
in the original basis. 
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Fig. 4.4(i). Sequences of F moves, braiding, and loop eliminations may be performed 
until the diagram has been reduced once more to a fusion tree without loops on a lattice 
of n sites. 

Much as the state of an anyonic system may be represented by a vector c", anyonic 
operators may be represented by a matrix M^^^ . Each value of a corresponds to a pair 
{aj,/ij} where is a possible charge of the central edge of the operator diagram [e.g. 03 
in Fig. 4.4(ii)], and /ij is a value of the degeneracy index associated with charge a^. We 
will denote the degeneracy of a, by Similarly, values of {3 correspond to pairs {oj, /ij} 
where Oj has degeneracy For degeneracy-preserving anyonic operators the charge 
indices and aj necessarily take on the same range of values, and z/q- = z/^^ when = aj. 
The values of may equivalently be calculated from either the fusion tree making up 
the top half or the splitting tree making up the bottom half of the operator diagram. 

A well-defined anyonic operator M must respect the (quantum) symmetry group of the 
anyon model, and consequently all entries in M^^ for which aj 7^ will be zero. However, 
in contrast with we do not require that = a^- = 1. When M is a degeneracy-preserving 
operator, matrix is therefore a square matrix of side length 



which may be organised to exhibit a structure which is block diagonal in the charge 
indices and a^, and for which the blocks are also square. As an example consider 
Fig. 4.5, which shows an operator acting on four Fibonacci anyons. An example matrix 
M^^ for an operator of this form is given in Table 4.1, from which the entries of Matcde 
can be reconstructed, e.g. Mririr = 3. 

4.3.2 Degeneracy- Changing Anyonic Operators 

We now introduce the second class of anyonic operators, which map states in some Hilbert 
space H into some other Hilbert space "H'. These operators may reduce or increase the 
degeneracy of any charge present in the spaces on which they act, and may even project 
out entire charge sectors by setting their degeneracy to zero. When these operators 




(4.5) 
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Figure 4.5 - An operator acting on four Fibonacci anyons. The values of the coefficients 
Mahcde may be specified as a block-diagonal matrix mJ^ , for example as in Table 4.1. 
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Table 4.1 - Matrix representation mJ^ for an example operator of the form shown in Fig. 4.5. 
Multi-index a corresponds to index pair {aj, /ij} and multi-index (3 corresponds to pair {aj, ^j}. 
Subject to an appropriate ordering convention for /Xj and fij, these indices may be related to 
the fusion tree labels a,b,c,d,e of Fig. 4.5 as shown. Note that as c is the charge on the 
central leg of Fig. 4.5, all nonzero entries of Ma^ satisfy Ui = aj = c. 
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are written in the conventional notation of Fig. 4.4, the fusion and sphtting trees will 
not be identical. Further, we may choose to allow combinations of degeneracies which 
do not naturally admit complete decomposition into individual anyons. For example, a 
degeneracy- changing operator may map a state on five Fibonacci anyons (having total 
degeneracies ui = 3, Ur = 5) into a state having degeneracies z/i = 2, Ur = 2. As these 
degeneracies do not admit decomposition into an integer number of nondegenerate anyons, 
it is necessary to associate an index Ui with the single open leg of the fusion tree. This 
index behaves identically to the vertex indices Ui of Fig. 4.1, serving to enumerate the 
different copies of each individual charge, and as with the vertex indices of Fig. 4.1, it is 
absorbed into the degeneracy index /Xj. 

As a further example, a state having degeneracies i/i = 4, = 4 could be associated with 
a fusion tree having either one leg, or two legs each with degeneracies ui = 0, Ut- = 2. 
Again, indices Ui would have to be associated with each open leg. 

Matrix representations of degeneracy-changing anyonic operators may also be constructed, 
and when they are written in block diagonal form, the matrices and their blocks may 
be rectangular rather than square. Degeneracy-changing anyonic operators therefore 
represent a generalisation of the degeneracy-preserving anyonic operators discussed in 
Sec. 4.3.1. It is worth noting that the presence of indices on the open legs of the 
fusion or splitting trees of an operator do not automatically imply that it is a degeneracy- 
changing anyonic operator: The defining characteristic of a degeneracy-preserving anyonic 
operator is that it maps a state in a Hilbert space H into a state in the same Hilbert space 
H, and consequently both the matrix as a whole and all of its blocks are square. Thus a 
degeneracy-preserving anyonic operator may act on states having additional indices Ui on 
their open legs, and the resulting state may be expressed in the form of the same fusion 
tree, with the same additional indices on the open legs. 

Operators which change degeneracies may represent physical processes which change the 
accessible Hilbert space of a system. As we will see in Sec. 4.5.1, they may also be used 
in tensor network algorithms as part of an efficient representation of particular states or 
subspaces of a Hilbert space, for example the ground state or the low energy sector of a 
local Hamiltonian. 
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Figure 4.6 - Diagrammatic representation of operators P^^^ (4-6) and P^^-' (4.7). The charges 
on all leaves are non-degenerate. 



This distinction between degeneracy- changing and degeneracy-preserving anyonic opera- 
tors is clearly seen with a simple example. Let {ip) be a state on six Fibonacci anyons. 
This state can be parameterised by a vector c"", which has five components. We now 
define two projection operators, P*-^^ and P*-^-* (Fig. 4.6), each of which acts on the fusion 
space of anyons ti and T2. Operator P^^^ is degeneracy-preserving, and projects into 
the subspace in which anyons ri and T2 fuse to the identity. Its matrix representation is 

where the first value of each multi-index corresponds to a charge of 1, and the second to 
a charge of r. Operator P^^^ performs the same projection, but is degeneracy- changing. 
Its matrix representation is written 

P(^V = (10). (4.7) 

Both operators perform equivalent projections, in the sense that 

(^|P«tp(i)|^) = (^|p(2)tp(2)|^). (4.8) 

When P(i) acts on lip) it leaves the Hilbert space unchanged, and hence the vector c'" 
describing state = P^^'^lip) is once again a five-component vector, although in an 
appropriate basis some components will now necessarily be zero. In contrast P*^^^ explicitly 
reduces the dimension of the Hilbert space, and the vector c"" describing state {tp") = 
P^'^^ip) is of length two, describing a fusion tree on only four Fibonacci anyons (as both ri 
and T2 have been eliminated). One consequence of this distinction is that while (P'^^-')^ = 
P^^\ the value of (P^^))^ is undefined. 
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4.4 Anyonic Tensor Networks 
4.4.1 Diagrammatic Notation 

The diagrammatic notation conventionally employed in the study of anyonic systems, and 
used here in Figs. 4.1 and 4.4, is well suited to the complete description of anyonic systems, 
as it provides a physically meaningful depiction of the entire Hilbert space. However, the 
number of parameters required for such a description grows exponentially in the system 
size, and because it is necessary to explicitly assign every index to a specific charge or 
degeneracy, specification of a tensor network rapidly becomes inconveniently verbose [for 
example see Fig. 4.4(iii)]. 

In the preceding Sections, we developed techniques whereby anyonic states and opera- 
tors could be represented as vectors and matrices, bearing only one or two multi-indices 
apiece. We now introduce the graphical notation which complements this description, 
and in which we will formulate anyonic tensor networks. Figure 4.7(i) gives the graphical 
representations of a state associated with a vector c", and of an operator M asso- 
ciated with a matrix M^^ . The circle marked c corresponds to the vector c", and the 
circle marked M corresponds to the matrix M^^. In general, grey circles correspond to 
tensors, and the number of legs on the circle corresponds to the number of multi-indices 
on the associated tensor. Each multi-index is also associated with a fusion or splitting 
tree structure, which is specified graphically. For reasons to be discussed shortly, we will 
require that no tensor ever have more than three multi-indices. As the legs of the grey 
shapes are each associated with a multi-index, they carry both degeneracy and charge 
indices. Consequently it is not necessary to explicitly assign labels to the fusion/splitting 
trees, as these labellings are contained implicitly in the degeneracy index (for example 
see Table 4.1, where specifying the values of {ai,/ij} and {aj,/ij} is equivalent to fully 
labelling the fusion and splitting trees of Fig. 4.5). 

The fusion or splitting tree associated with a particular multi-index may be manipulated 
in the usual way by means of braids and F moves, recalling that each component of the 
tensor is associated with a particular labelling of the fusion and splitting trees via the 
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Figure 4.7 - (i) Diagrammatic representation of a state and two-site operator M expressed 
in terms of degeneracy indices, (ii) Application of M to state \ip). Grey shapes represent 
tensors with charge and degeneracy multi-indices, with each leg of the shape corresponding 
to one charge and degeneracy index pair. These diagrams represent the same state, operator, 
and process as Fig. 4.1(i) and Fig. 4.4(ii)-(iii). 
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corresponding values of the multi-indices. Manipulations performed upon a particular 
tree thus generate unitary matrices which act upon the multi-index that corresponds to 
the labellings of that particular tree. 

The application of an operator to a state is, unsurprisingly, performed by connecting 
the appropriate diagrams, as shown in Fig. 4.7(ii). For operators of the type discussed 
in Sec. 4.3.1, the outcome is necessarily a new state in the same Hilbert space, which 
consequently can be described by a new state vector c'°, as shown. However, in general 
an operator M will not act on the entire Hilbert space of the system, and so will be 
described by a tensor constructed on the fusion space of some subset of lattice sites, and 
not on the system as a whole. Operator M acting on state {ip) in Fig. 4.7(ii) is an example 
of this. Because c° describes a six-site system but is constructed on the fusion space 
of two sites, the multi- indices of M^,^ span a significantly smaller Hilbert space than that 
of c° and we cannot simply write 

c'^ = c^Mj (4.9) 

(using Einstein notation, where repeated multi-indices are assumed to be summed). In- 
stead, we must understand how to expand the matrix representation of an operator on 
some number of sites x, to obtain its matrix representation as an operator on x' sites, 
where x' > x. 

4.4.2 Site Expansion of Anyonic Operators 

The multiplicity tensor A*"^^ describes the fusion of two charges without degeneracies. It 
is easily extended to incorporate degeneracies of the charges, and we will denote this 
expanded multiplicity tensor N2p^ where multi-indices a, (3, and 7 are associated with 
the pairs {a,fia}, {6, /i^,}, and {c, /ic} respectively, and for given values of a, (3, and 7, 
u runs from 1 to A^^. The degeneracies associated with charges a, b, and c are denoted 
Ua, z/f,, and z/c respectively. As with /Iq, /if,, and fic, there is an implicit additional index 
on each degeneracy representing the edge of the tree on which charge x resides. The 
values of z/^ and Uh may be chosen arbitrarily (for example, i'a\a=i niay differ from z/b|b=i). 
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Table 4.2 - Construction of N^/Bu ^ fusion vertex for Fibonacci anyons. In this example 
a may take charges 1 and r each with degeneracy 1, and b may take charges 1 and r with 
degeneracies 1 and 2 respectively. By Eq. (4.2), charge c may therefore take values 1 and r 
with degeneracies 3 and 5 respectively. A correspondence between the values of multi-index 7 
and of multi-indices a and (3 is established in some systematic manner, with each assignation 
satisfying c £ a x b, and for Fibonacci anyons the index u is trivial as all multiplicities N^^^ 
are zero or one. An example assignation is shown in the table. The corresponding entries of 
N are then set to 1, with all other entries zero. For example, the fourth row indicates that 

^\l,l)(r,l)l - ^■ 

but the degeneracies associated with the values of c must satisfy 

= 5^ ^a^6iV„^, (4.10) 
a,b 

in accordance with Eq. (4.2). When this constraint is satisfied, every quadruplet of indices 
{a, fia,b, fib} corresponding to a unique pair of choices for a and (3 may be associated 
with N^fj distinct pairs of indices {c, /ic} for each c G a x 6. These pairs {c, fic} are 
enumerated by the additional index u. This defines a 1:1 mapping between sets of values 
on {a, ^a,b, ^b,u} and pairs {c, /ic}, and we set the corresponding entries in A^2/3« 
with all other entries being zero. A simple example is given in Table 4.2. 

By virtue of their derivation from A^^, the object AL^ and its conjugate iV^"'^" represent 
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application of the anyonic fusion rules, and may be associated with vertices of the splitting 
and fusion trees. Under the isotopy invariance convention there is an additional factor of 
[df./ {dadh)]^ associated with the fusion of charges a and h into c, where d^ is the quantum 
dimension of charge x, and similarly for splitting, but we will account for these factors 
separately (see Sec. 4.7.1). Thus constructed, the tensors N satisfy NI^^N'^'^^^ = S^'^. 

When used as a representation of the fusion rules, the generalised multiplicity tensor N2i3u 
and its conjugate iV^"''" permit us to increase or decrease the number of multi-indices on 
a tensor in a manner which is consistent with the fusion rules of the quantum symmetry 
group. This process is reversible provided the symmetry group is Abelian or, for a non- 
Abelian symmetry group, provided the total number of multi-indices on the tensor does 
not at any time exceed three. In constructing and manipulating a tensor network for a 
system of anyons, we will require only objects which respect the fusion rules of the anyon 
model. It is a defining property of such objects that when the number of multi-indices 
they possess is reduced to 1 by repeated application of and A^"'', non-zero entries may 
be found only in the vacuum sector. We imposed this requirement for states in Sec. 4.2, 
and it is equivalent to the restriction we imposed on anyonic operators in Sec. 4.3.1. In 
Singh et al. (2010a) (Sec. 3.1 of this Thesis), an equivalent condition was observed for 
tensors remaining unchanged under the action of a Lie group, and these tensors were 
termed invariant. When working with invariant tensors, we may separately evaluate the 
components of the tensors acting on the degeneracy spaces (e.g. the nonzero blocks of 
M^^), and the factors arising from loops and vertices of the associated spin network (see 
Sec. 4.7.1 for details). This property greatly simplifies the contraction of pairs of tensors. 

In addition to increasing or decreasing the number of legs of a tensor, we may also use A^ 
to "raise" the matrix representation of an operator from the space of x sites to the space 
of (x + x') sites. This is shown in Fig. 4.8, and the matrix representation of the raised 
operator is given by 

M'/ = M/iVtr iVL (4.11) 

where multi-index e describes the fusion space of all sites in {x + x') but not in x. Because 
the numeric factors associated with loops and vertices (and braiding where applicable) 
are handled separately, no factors of quantum dimensions appear in Eq. (4.11). 
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Figure 4.8 - "Raising" of an operator M from sites x to sites x + x': (i) Operator M defined 
only on sites denoted x. (ii) Resolutions of the identity are inserted above and below M, 
being constructed from tensors N and W . The central portion of this diagram is identified 
as corresponding to the new matrix Ma^ which describes M on x + x'. (iii) Loop and vertex 
factors in the central region are evaluated separately and eliminated, (iv) The tensor network 
corresponding to the new central portion is contracted. The N and W tensors outside the 
central region become vertices of the fusion and splitting trees associated with Ma^ . Together 
the trees and the matrix ML^ constitute the raised version of M. 
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To act an operator M on a state I?/') in the matrix representation, we therefore connect 
the diagrams for M and lip), ehminate all loops, and then raise the matrix representation 
of the operator M using Eq. (4.11), repeatedly if necessary, until the resulting matrix 
M'q,^ may be applied directly to the state vector c". Similarly it is possible to combine 
the matrix representations of operators, by connecting their diagrams appropriately, elim- 
inating loops, and performing any required raising so that both operators act on the same 
fusion space. Their matrix representations can then be combined to yield the matrix 
representation of the new operator: 

^(1X2) (S ^ ^(1) ,^(2) (S ^4_^2) 

and the fusion/splitting tree associated with this new operator is obtained as shown in 
Fig. 4.8. 

Note that as yet, we have not described how two objects may be combined if their multi- 
indices are both up or both down, and are connected by a curved line. To contract such 
objects together, it is necessary to understand how bends act on the central matrix of 
an operator. Once this is understood, the bend can be absorbed into one of the central 
matrices, so that the connection is once again between an upper and a lower multi-index 
as in Eq. (4.12). This process is described in Sec. 4.4.3. 

4.4.3 Manipulation of Anyonic Operators 

As observed in Sec. 4.4.2, when we describe a system entirely in terms of objects invariant 
under the action of the symmetry group, we may account separately for the numerical 
normalisation factors associated with the spin network. However, as well as affecting these 
numerical factors, transformations of the fusion or splitting tree of an anyonic operator 
will typically also generate unitary matrices which act on the matrix representation of 
the operator. These matrices respect the symmetry of the anyon model, and thus can 
be written as block-diagonal matrices where each block is a unitary matrix acting on a 
particular charge sector. In terms of the diagrammatic notation of Sec. 4.4.1, F moves and 
braids therefore result in the insertion of a unitary matrix, as shown in Fig. 4.9. These 
matrices, whose entries are derived from the tensors {Fl}^'^)(^euv){fu'v') and Rf' respectively. 
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Figure 4.9 - (i) F move, and (ii) braiding, performed on a section of fusion tree in the 
diagrammatic notation of Sec. 4.4.1. 

are raised if required, as described in Sec. 4.4.2, and then contracted with M^^, the matrix 
representation of the operator. To compute the unitary matrices involved, it suffices to 
recognise that F moves and braids are unitary transformations in the space of labelled 
tree diagrams. Identifying the leg on which the unitary matrix is to be inserted, the 
relevant region of the space of labelled diagrams is then enumerated by the multi-index 
which can be associated with this leg (compare Fig. 4.2). 

Braiding is of particular importance when working in two dimensions, as an operator 
will necessarily be defined with respect to some arbitrary linear ordering of its legs, and 
when manipulating a tensor network it may be necessary to map between this original 
definition and other equivalent definitions, corresponding to different leg orderings. For 
example, let M be a four-site anyonic operator as shown in Fig. 4.10(i), which we wish 
to apply to a 2D lattice. For the indicated linearisation of this lattice, application of M 
will require braiding as shown in Fig. 4.10(ii). By evaluating the unitary transformations 
corresponding to these braids and absorbing them into M^,^, we may define a new operator 
M' which acts directly on the linearised lattice without any intervening manipulations of 
the fusion/splitting trees. 

We will also frequently wish to deal with tensor legs which bend vertically through 180°. 
If working with an anyon model that has non-trivial Frobenius-Schur indicators, then 
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(i) 



12 3 4 



(ii) 



M 



M 




1 2 3 




(iii) 



13 2 4 




Figure 4.10 - (i) An operator M acting on sites on a 2D lattice is defined with respect to 
some arbitrary linear ordering of these sites, (ii) When manipulating the tensor network, it may 
on occasion be computationally convenient for the lattice to be linearised according to some 
alternative linearisation scheme. In this example, the imposed linearisation scheme is indicated 
by the dotted line. To apply M to a different linearisation of the lattice may require braiding. 
The orientation of the braids can be determined by putting the fusion tree of (i) onto the 2D 
lattice, then smoothly deforming the lattice into a chain in accordance with the linearisation 
prescription, (iii) The unitary matrices corresponding to the required F moves and braiding 
operations may be absorbed into M, defining a new operator M' on the linearised lattice. 
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indicator flags must be applied to all bends. Like F moves and braiding, the reversal of a 
Frobenius-Schur indicator flag is a unitary transformation, and once again this leads to 
the introduction of a unitary matrix which can be absorbed into a nearby existing tensor. 
However, we may wish to perform other operations on bends, such as absorbing them into 
fusion vertices or the central matrices of anyonic operators. We may also need to move 
a matrix M^,^ across a bend. We must therefore develop the description of bends in the 
new diagrammatic formalism. 

In Bonderson (2007) a prescription for absorbing bends into fusion vertices is given in 
terms of tensors {Af')^^ and {B^^)^^, derived from the F moves, and corresponding to 
clockwise and counter-clockwise bends respectively. The absorption of a clockwise or 
counterclockwise bend into a fusion vertex is reproduced in Fig. 4.11(i), and results in a 
vertex fusing upward- and downward-going legs. We now assign new tensors (iV^^)^"^ 
and (iV*-"-^^)^^" to such vertices, such that writing these transformations in the notation 
of Sec. 4.4.1 is trivial. This is shown in Fig. 4.11(ii). 

Explicit expressions for the new vertex tensors (A^*"^)"'^ and (A^'-'*-^^)"'' may be obtained 
by recognising that Fig. 4.11(i) describes the action of unitary transformations on AT^"^**. 
When the bend is counterclockwise, the corresponding unitary matrix is derived from 
(y4j?*)„t,, and when the bend is clockwise, the unitary matrix is derived from {Bf')uv We 
will denote these unitary matrices A^^ and B^^ respectively. We then have 

(iVCW)t- = A^'N^r^ep (4.13) 
(^ccw)t/3n ^ i?^^ivtf"5,,. (4.14) 

and conjugation describes equivalent vertices and N'~^^'^ when a bend is absorbed 

into a splitting tree. 

Knowing how the absorption of bends acts on a vertex tensor, we may readily infer how 
the same process acts on the matrix representation of an operator. In Fig. 4.11(iii) we see 
a bend absorbed into the matrix M^^, resulting in a new object with two lower multi- 
indices, M'^p. First we exploit the freedom to introduce fusion with the trivial charge 
(denoted I), with degeneracy 1. The corresponding A^''" object takes only one value on its 
upper left multi-index, and is fully defined by iV^^^^ = 5^ . Absorbing the bend into this 
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Figure 4.11 - Vertical bending of legs (i) in the standard diagrammatic notation, and (ii) in 
the diagrammatic notation of Sec. 4.4.1. White triangles represent Frobenius-Schur indicator 
flags, (iii) Legs on the matrix representations of states and operators may also absorb bends. 
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(i) ^ (ii) 

u 

Figure 4.12 - Opposing pairs of Frobenius-Schur indicators (i) on a pair of bends equivalent 
to the identity, and (ii) on a pair of bends such as might be used when computing a quantum 
trace, (iii) As an anyon model can always be specified such that the Frobenius-Schur indicators 
are ±1, reversing a pair of contiguous opposed Frobenius-Schur indicator flags is always free. 




Figure 4.13 - Moving a matrix across a bend in a tensor network diagram, (i) Initial diagram, 
(ii) Bends are absorbed into the matrix, (iii) New bends are introduced, in accordance with 
Fig. 4.12(i). (iv) A pair of contiguous, opposed Frobenius-Schur indicators are reversed, as per 
Fig. 4.12(iii). The initial and final matrices M and M' are related as specified in Eq. (4.16). 



fusion vertex as per Eq. (4.13) yields A^^6g'6,^ = A^p, which may be then combined with 
Mj to give 

M^^ = M^-'A,^. (4.15) 

In conjunction with the relationships given in Fig. 4.12, this gives us the ability to move 
a matrix past a bend. An example of this is given in Fig. 4.13, for which M and M' are 
related according to 

^ = A„^M,^x/5t^^ (4.16) 

where >i/ represents reversal of the Frobenius-Schur indicator flag on the lower bend. 
Finally, bending may also allow more efficient contraction of pairs of anyonic operators, 
as shown in Fig. 4.14. 

Having described the action of bends, it is customary also to introduce a second type of F 
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(i) 



V (ii) 



(ii) 




(i) 




(ii) A 



(ii) 
..>... 











(b) 










Figure 4.14 - The use of bends may permit the more efficient contraction of pairs of anyonic 
operators. In the sequence of events marked (i), operator A is first raised to the space of 
three sites then contracted with B. In sequence (ii) the operators are instead contracted 
using bends. For many anyon models the latter approach offers a significant computational 
advantage. 
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a, fl. 



(a^,Ui,U2)(af,u^,u^) 



a, a. 



Figure 4.15 - Now that we may bend legs up and down it is customary to introduce a further 
type of F move, derived by applying bends to the one presented in Fig. 4.1(iii). 



move which is described by the tensor (F^^^^^) {a5UiU2)iaeu3U4) (Fig. 4.15). This tensor may be 
derived from (i^a4'''"'')(a5«i«2)Kn3«4) by bending, and as with (i^a4'''"")(a5«i«2)(a6«3«4) these 
F moves perform a transformation of the fusion tree, accompanied by the introduction of 
a unitary matrix which can be absorbed into the matrix representation of the operator. 
These unitary matrices correspond to the consecutive apphcation of a bend, an F move 
of the original type, and a second bend whose action is the inverse of the first. 

4.4.4 Constructing a Tensor Network 

Now that we have developed a formalism for anyonic tensors, we may convert an exist- 
ing tensor network algorithm for use with anyons. First, the tensor network must be 
drawn in such a manner that every leg has a discernible vertical orientation. Although 
these orientations may be changed during manipulation of the tensor network, an initial 
assignment of upward or downward direction is required. Second, all tensors must be 
represented by entirely convex shapes, such as circles or regular polygons. For existing 
tensor network algorithms such as MERA and PEPS, this requirement is trivial. However, 
it is conceivable that future algorithms might involve superoperator-type objects whose 
graphical representations interleave upward- and downward-pointing legs. Concavities on 
these objects may be eliminated by replacing some of their upward-pointing legs with 
downward-pointing legs (or vice versa), followed by a bend [Fig. 4.16(i)-(ii)]. A similar 
treatment may be applied to any superoperators which arise during manipulations of the 
tensor network, introducing a pair of bends as in Fig. 4.12(i) and then absorbing one into 
the matrix representation of the object. 

If working with an anyon model that has non-trivial Frobenius-Schur indicators, then 
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(i) 



M 



(ii) 



> 



n 



< 



(iii) 



> 



M 




Figure 4.16 - Construction of an anyonic tensor corresponding to a normal tensor with more 
than three legs, (i) The original tensor, (ii) If required, any concavities are eliminated by 
introducing bends, (iii) Frobenius-Schur indicators are assigned to the bends, (iv) Directions 
are assigned to all legs, consistent with the rest of the network, (v) Legs are collected together 
into fusion and splitting trees. The central object, representing degrees of freedom of the 
tensor, now has less than four legs, (vi) If desired, bends can be re-absorbed into the fusion 
and splitting trees. 



indicator flags must be applied to all bends. Initial choices are a matter of convenience, 
and it is frequently possible to assign these indicators in opposed pairs, as shown in 
Fig. 4.12. If these paired indicators are not flipped or are only flipped in adjacent opposed 
pairs during subsequent manipulations of the tensor network, then they may frequently 
be left implicit. 

Next, if there exist charges in the anyon model which are not self-dual, a direction (rep- 
resented by a solid arrow) must be assigned to every multi-index. Any tensor with more 
than three legs (e.g. M in Fig. 4.16) is then replaced by a trivalent tensor network con- 
sisting of a core object, e.g. M^,^, which contains the free parameters of the tensor, and 
as many copies of or N'^ as are required to provide the correct output legs. These 
tensors A^, A^^ correspond to vertices in the fusion and splitting trees associated with 
M^'^, yielding the corresponding anyonic tensor. Objects with three legs or less can be 
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directly identified with an anyonic tensor object carrying the appropriate number of in- 
dices (i.e. three multi- indices and a vertex index u), though for consistency with the 
methods described in Sections 4.4.2 and 4.4.3 we point out that it is possible to similarly 
replace three-legged objects with anyonic operators consisting of a central matrix M^^ 
and a fusion or splitting vertex, if desired. 

Any bends introduced earlier may now be reabsorbed, so that some vertices now cor- 
respond to (A^W)^ (A^cw)^ (A^w^t^ and (A^^w^t. This step, however, is optional as 
it may be more convenient for subsequent manipulations of the tensor network if the 
bends are left explicit. The anyonic tensors are then connected precisely as in the original 
Ansatz. 

Manipulations of the anyonic tensor network are equivalent to those performed on the 
spin version of the Ansatz, differing only in that the degrees of freedom of the tensor 
network are now expressed entirely by the at-most-trivalent central objects, and certain 
topological elements such as braids and vertical bends must be accounted for in accordance 
with the prescriptions of Sec. 4.4.3. These changes may naturally imply minor changes to 
the manipulation algorithms, and we will see examples of this in the ID MERA. Similar 
considerations will apply to other tensor network algorithms. 

Our construction of an anyonic tensor network draws upon two important elements which 
have previously been observed in other, simpler, physical systems: 

1. Tensors in the Ansatz exhibit a global symmetry, which may be non-Abelian. Ex- 
ploiting a non-Abelian symmetry requires that the Ansatz be written in the form 
of a trivalent tensor network. This has previously been observed and implemented 
for non-Abehan Lie group symmetries such as SU(2) (Singh et al., 2010a, Sec. 3.1 
of this Thesis; Singh and Vidal, in preparation). 

2. Tensors in the Ansatz must be able to account for non-trivial exchange statistics. 
This has previously been observed in the simulation of systems of fermions (Barthel 
et al., 2009; Corboz et al., 2010a,b; Gu et al., 2010; Kraus et al., 2010; Pineda et al., 
2010; Pizorn and Verstraete, 2010; Shi et al., 2009), where efficient implementation 
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of particle statistics can be achieved through the use of "swap gates" (Barthel et al., 
2009; Corboz et ah, 2010b; Pineda et al., 2010; Pizorn and Verstraete, 2010). 

In both cases, anyonic tensor networks extend the concepts introduced in previous work. 
The symmetry structure of an anyon model may be a quantum group, for example a mem- 
ber of the series S'f/(2)fc, k G Z+, rather than having to be a Lie group, and this permits 
representation of non-Abelian anyonic systems whose Hilbert space does not admit decom- 
position into a tensor product of local Hilbert spaces. Similarly, anyonic braiding may be 
implemented using a generalisation of the fermionic "swap gate" formalism. When braid- 
ing, particle exchange may introduce transformation by a unitary matrix rather than by 
a sign, and efficient implementation of the resulting swap gates is particularly important 
for the simulation of 2D systems. 

Although anyonic systems pose a number of unique challenges, we see that these are 
addressed by developments based on existing techniques, and we therefore anticipate that 
the resulting generalisations of existing tensor network Ansatze should still be capable of 
accurately representing the states of an anyonic system. 

4.4.5 Contraction of Anyonic Tensor Networks 

The techniques described in Sees. 4.4.2 and 4.4.3 (F moves, braids, bending of legs, 
elimination of loops, diagrammatic isotopy, flipping of Frobenius-Schur indicator flags, 
and the use of N^'^^ tensors) suffice to contract any network of anyonic tensors written 
in the form of matrices with degeneracy indices, and unlabelled trees. Through careful 
application of these techniques, and avoiding at all times processes which would yield a 
tensor with more than three legs, the matrix representations of any pair of contiguous 
tensors in a network may always be brought into conjunction such that their multi-indices 
can be contracted in the manner of Eq. (4.12), and any tensor network may be contracted 
by means of a sequence of such pairwise contractions. 

That a tensor network may represent a system of anyons in this way is possible because 
throughout the anyonic tensor network, each value of a degeneracy index is associated 
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with a specific labelling of the corresponding unlabelled tree. Consequently it is always 
possible to fully reconstruct any operation in terms of the more verbose representation of 
Fig. 4.4. 

An anyonic tensor network is therefore fully specified merely by the unlabelled tree (with 
Frobenius-Schur indicator fiags if required), and the values and locations of the matrix 
representations of its tensors, written in the degeneracy index form. 

4.5 Example: The ID MERA 
4.5.1 Construction 

To construct an anyonic MERA for a ID lattice with n sites, where n satisfies n = 
2 X 3^, k G Z+, we begin with a "top" tensor on a two-site lattice Cr whose matrix 
representation is of a computationally convenient size. [The top tensor is named for its 
position in the usual diagrammatic representation of the MERA, where diagrams with 
open legs at the bottom correspond to a ket. For anyons the converse convention applies, 
and consequently in Fig. 4.17(i) the "top" tensor is ironically located at the bottom.] 

To each leg of the top tensor, we now append an isometry [Fig. 4.17(ii)]. The matrix 
representations of the isometrics consist of rectangular blocks, as described in Sec. 4.3.2, 
and we choose isometrics whose fusion trees have three legs, so as to construct a ternary 
MERA (Evenbly and Vidal, 2009a). Next, disentanglers are applied above the isometrics. 
For periodic boundary conditions this must be performed in a manner which respects the 
anyonic braiding rules, as shown in Fig. 4.17(iii). We identify the open legs of the resulting 
network as the sites of a lattice Cr-i, and the rows of disentanglers and isometrics may be 
understood as a coarse-graining transformation taking a finer-grained lattice Cr-i into a 
coarser-grained lattice Cr, similar to the standard MERA. Note that the geometry of the 
periodic lattice is refiected by the connections of the disentanglers. Specifically, whether 
the outside legs are braided over or under the other lattice sites refiects whether the lattice 
closes towards or away from the observer. 
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Figure 4.17 - Construction of a ID ternary MERA on a periodic lattice from anyonic operators, 
(i) The "top" tensor, T. (ii) Isometries, w. (iii) Disentanglers, u. The fusion tree representing 
an anyonic state (or ket) is usually drawn with the lattice sites at the top, so this MERA 
has been constructed "upside down" when compared with the diagrams in Evenbly and Vidal 
(2009a) and Vidal (2010). This is unimportant, and we could equally well have decided to 
follow the convention usually adopted in tensor network algorithms, labelled the tensors in (i)- 
(iii) by T"^ , , and u"^ , and identified diagram (i)-(iii) as a bra. (iv) Structure of a 2-site 
term in the Hamiltonian, h, or a 2-site reduced density matrix, p. (v) Disentanglers and (vi) 
isometries satisfy the relationships it'^u = I, w'^w = I. 
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The application of anyonic isometries and disentanglers is now repeated k times [Fig. 4.17(i)- 
(iii) corresponds to = 1], until the Ansatz has n legs. The final row of isometries should 
be chosen such that each of their upper legs have the same charges and degeneracies as the 
sites of the physical lattice Cq, and the open legs above the last row of disentanglers are 
identified with the physical lattice. For coarse-grained lattices Ci to Cr, the dimensions 
of the lattice sites correspond to the lower legs of the isometries and are chosen for com- 
putational convenience, subject to the requirement that each charge sector is sufficiently 
large to adequately reproduce the physics of the low-energy portion of the Hilbert space. 
For all other legs, their charges and degeneracies are determined by requiring consistency 
with Eq. (4.2). Initial choices of which charges to represent on the "top" tensor and on 
the lower legs of the isometries, and with what degeneracies, must be guided either by 
prior knowledge about the physical system, or by balancing computational convenience 
against the inclusion of a broad and representative range of possible charges. When used 
in a numerical optimisation algorithm, the choice of relative weightings for the different 
charge sectors may often be refined by examination of the spectra of the reduced density 
matrices on the coarse-grained lattices, after initial optimisation of the tensor network is 
complete. 

This concludes construction of the MERA for a state on a finite, periodic ID anyonic 
lattice. That this tensor network does represent an anyonic state is easily seen by se- 
quentially raising tensors, performing F moves, and combining tensors, until the entire 
network is reduced to a single vector whose length is equal to the dimension of the phys- 
ical Hilbert space, and an associated fusion tree. These then represent the state of the 
system as per Eq. (4.4). The structure of this tensor network closely resembles that of 
the normal MERA, according to the identifications given in Fig. 4.17, and consequently 
we anticipate that it will share many of the same properties, including the ability to re- 
produce polynomially decaying correlators in strongly correlated physical systems. Open 
lattices may also be easily represented by omitting the braided disentanglers at the edge 
of the diagram. 

We also note that in common with the MERA for spins, the anyonic MERA may be 
understood as a quantum circuit, although one which carries anyonic charges in its wires. 
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Any junction in the fusion/sphtting trees may be associated with a or A^^ tensor, and 
the entire network may be considered as the apphcation of a series of gates to a Hilbert 
space of fixed dimension beginning mostly (or entirely, if the top tensor is considered to 
be the first gate) in the vacuum state, with individual gates introducing entanglement 
across some limited number of wires. 

4.5.2 Energy Minimisation 

The anyonic MERA can be used as a variational Ansatz to compute the ground state 
of a local Hamiltonian. The Hamiltonian is introduced as a sum over nearest neighbour 
interactions, each term having the form of Fig. 4.17(iv), and optimisation of the tensor 
network is carried out in the usual manner (Evenbly and Vidal, 2009a). Also as per 
usual, Hamiltonians involving larger interactions, such as next-to-nearest neighbour, can 
be accommodated by means of an initial exact n-into-one coarse- graining of the physical 
lattice. 

As in Evenbly and Vidal (2009a), optimisation of the MERA then consists of repeatedly 
lifting the Hamiltonian from Cq to the coarse-grained lattices, updating their isometrics 
and disentanglers, and lowering the reduced density matrix, or the top tensor and its 
conjugate. When lifting the Hamiltonian or lowering the reduced density matrix, then 
the diagrams in Evenbly and Vidal (2009a) taken in conjunction with the key given in 
Fig. 4.17 serve to describe networks of anyonic operators which, when contracted to a 
single operator, yield the lifted form of the Hamiltonian or lowered form of the reduced 
density matrix respectively. Similarly, when optimising disentanglers or isometrics, the 
diagrams of Evenbly and Vidal (2009a) and the identifications in Fig. 4.17 indicate how 
to construct an anyonic operator which constitutes the environment of the anyonic oper- 
ator being optimised. However, once the admissible ranges of charges and degeneracies 
on each leg have been fixed, the only optimisable content of an anyonic operator is its 
matrix representation. Consequently, the fusion and splitting tree contributions should 
be evaluated and absorbed into the operator and its environment, reducing them both to 
their matrix representations, denoted M and E respectively (see Fig. 4.18). If the sin- 
gular value decomposition of E is written E = USW"^ , then the updated matrix content 
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(i) 




(ii) 



■■>■ 




(iii) 



(iv) 



> 



> 



Figure 4.18 - (i) Anyonic operator E constitutes the environment of operator M. Factors 
arising from the fusion and splitting trees should be evaluated and absorbed into matrices E 
and M, following which (ii) matrix E constitutes the environment of matrix M. After (iii) 
updating the matrix M to M' , (iv) the fusion and splitting trees of M should be reinstated, 
the numerical factors associated with this process being the inverse of the fusion tree factors 
previously absorbed into matrix M. Frobenius-Schur flags in (i)-(ii) are represented by white 
triangles, and are not to be confused with the black arrows which indicate the orientation of 
lines in the fusion/splitting trees. 



M of the anyonic operator being optimised is given by —WU\ minimising the value of 
Ti{EM) subject to the usual constraint for disentanglers and isometries that MM"!" = I 
[Fig. 4.17(v)-(vi)]. The fusion/splitting tree content of the operator can then be restored, 
along with any appropriate numerical factors that may be required. 

As with the standard MERA, the "top" tensor is constructed by diagonalising the total 
Hamiltonian on the most coarse-grained lattice, f^tot on ^^t- As £^ is a two-site lattice, the 
total Hamiltonian f^tot is a sum of two terms, H12 and H2i- For the translation- invariant 
anyonic MERA, we may formally define H21 in terms of H12 as shown in Fig. 4.19, and 
the top tensor T (together with any factors arising from the chosen normalisation scheme) 
then corresponds to the lowest-energy eigenstate of Htot- 
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Figure 4.19 - Definition of H21 in terms of H12, on the most coarse-grained lattice {Cr) of 
the translation-invariant periodic ternary MERA. Lattice Cr is a two-site periodic lattice. 



4.5.3 Scale- Invariant MERA 



Having identified the anyonic counterparts of the tensors of the standard MERA, and 
described how these tensors may be hfted, lowered, and optimised, the algorithm for 
the scale-invariant MERA described in Pfeifer et al. (2009) may also be implemented 
for anyonic systems, simply by applying the dictionary of Fig. 4.17 and the techniques 
described in Sec. 4.5.2. As with optimisation of u and w, the computation of the top 
reduced density matrix (which is a descending eigenoperator of the scaling superoperator 
with eigenvalue 1) may be understood as a calculation of the matrix component pj^ of 
the reduced density matrix p. The ascending eigenoperators of the scaling superoperator, 
or local scaling operators of the theory, may also be computed in this manner. 
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4.5.4 Results 



To demonstrate the effectiveness of the anyonic generahsation of the MERA, we apphed it 
to a ID critical system of anyons whose physical properties are already well known: The 
golden chain (Feiguin et al., 2007). This model consists of a string of Fibonacci anyons 
subject to a local interaction. Fibonacci anyons have only two charges, 1 (the vacuum) 
and r, and one non-trivial fusion rule (r x r — )■ 1 + r). The simplest local interactions for 
a chain of Fibonacci r anyons are nearest neighbour interactions favouring fusion of pairs 
into either the 1 channel (termed antiferromagnetic, or AFM), or the r channel (termed 
ferromagnetic, or FM). Both choices correspond to critical Hamiltonians, associated with 
the conformal field theories A^(4, 3) and A^(5, 4) for AFM and FM couplings respectively. 
Individual lattice sites are each associated with a charge of r. 

The AFM and FM Hamiltonians act on pairs of adjacent Fibonacci anyons. On a pair of 
lattice sites each carrying a charge of r, the matrix representations of the AFM and FM 
Hamiltonians are written 



where a multi-index value of 1 corresponds to the vacuum charge, 2 corresponds to r, 
and the charges are non-degenerate. We optimised a scale-invariant MERA on the golden 
chain for each of these Hamiltonians, and computed local scaling operators using the 
tensor network given in Fig. 4.20. The operators calculated using this diagram may be 
classified according to the values of the charge labels yi and y2-, and the scaling dimensions 
and conformal spins which we obtained are given in Tables 4.3 and 4.4, and Fig. 4.21. 

Comparison of the AFM case with existing results in the literature show that the scaling 
dimensions obtained when yi = y2 correspond to those obtained when studying a system 
of anyons with a toroidal fusion diagram (Feiguin et al., 2007). For a system of anyons 
on the torus it is possible to define an additional topological symmetry (Feiguin et al., 
2007) and classify local scaling operators according to whether or not they respect this 
symmetry. Operators satisfying yi = 1/2 = 1 correspond to those which respect the 
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Table 4.3 - Scaling dimensions for Fibonacci anyons with antiferromagnetic nearest neigbour 
interactions on an infinite chain. Numerical values were computed using an anyonic MERA 
with maximum degeneracies for charges 1 and r of 3 and 5 respectively (denoted x = [3,5]), 
and are grouped according to their classification by the values of yi and y2 in Fig. 4.20. 
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Table 4.4 - Scaling dimensions for Fibonacci anyons with ferromagnetic nearest neigbour 
interactions on an infinite chain. Numerical values were computed using an anyonic MERA 
with maximum degeneracies for charges 1 and r of 3 and 5 respectively (denoted x = [3,5]), 
and are grouped according to their classification by the values of yi and y2 in Fig. 4.20. 
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■' ' ' '. 




Figure 4.20 - Determination of eigenoperators {(/)) and associated scaling dimensions (A) for 
the one-site scaling superoperator of the anyonic ID MERA. Eigenoperators may be classified 
according to the charges on edges yi and y2- One interpretation of these labels is that, in 
addition to the sites of the ID lattice, there may exist free charges lying in front of and behind 
the anyon chain. The labels yi and ?/2 then represent the transfer of charge between these 
regions and the ID lattice. 

topological symmetry, and those satisfying yi = y2 = t do not. We will discuss the 
interpretation of the different sectors and their relationship to anyons on the torus in a 
forthcoming paper (Pfeifer and Vidal, 2010). 

When yi ^ y2 the scaling operators obtained are chiral, with those obtained from y^ = 
1, 2/2 = T and yi = T, ?/2 = 1 believed to form conjugate pairs. 

4.6 Summary 

Numerical study of systems of interacting anyons is difficult due to their non-trivial ex- 
change statistics. To date, study of these systems has been restricted to exact diagonalisa- 
tion. Matrix Product States (MPS) for ID systems, or special-case mappings to equivalent 
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Figure 4.21 - Scaling dimensions of leading primary operators and their descendants, com- 
puted for (i) antiferromagnetic and (ii) ferromagnetic local Hamiltonians (4.17) on the golden 
chain. Results are grouped into conformal towers, with a slight horizontal spread introduced 
to show the degeneracies of the descendant fields. A circled cross indicates a primary field, 
and a plain cross indicates a descendant. Dashed lines indicate values predicted from CFT. 
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spin chains. This paper shows how any tensor network Ansatz may be translated into a 
form apphcable to systems of anyons, opening the door for the study of large systems of 
interacting anyons in both one and two dimensions. As an example, this paper demon- 
strates how the MERA may be implemented for a ID anyonic system. This Ansatz is 
particularly important as many ID systems of anyons are known which exhibit extended 
critical phases (see e.g. Feiguin et al., 2007; Trcbst et al., 2008a,b). The structure of 
the MERA is known to be particularly well suited to reproducing long range correlations, 
and the scale-invariant MERA has the additional advantage of providing simple and direct 
means of computing the scaling dimensions and matrix representations of local scaling 
operators. 

We applied the scale invariant MERA to infinite chains of Fibonacci anyons under antifer- 
romagnetic and ferromagnetic nearest neighbour couplings, and identified a large number 
of local scaling operators. Our results for the scaling dimensions are in agreement with 
those previously obtained by exact diagonalisation of closely related systems, and for the 
relevant primary fields they are within 2.8% of the theoretical values obtained from con- 
formal field theory. We thus demonstrate that an anyonic MERA with x = [3, 5] permits 
conclusive identification of the relevant conformal field theory, and gives a level of accu- 
racy comparable to that of the scale invariant MERA on a spin chain (Table 2.1 of this 
Thesis; Pfeifer et al., 2009). 

The anyonic generalisation of the ID MERA presented here is useful in its own right, but 
the greatest significance of the approach described is that it is equally applicable to 2D 
tensor network Ansatze, and hence opens the door to studying the collective behaviour 
of large systems of anyons in two dimensions by numerical means, in situations where 
analytical solutions may not be possible. 

Note — Simultaneous with the work presented in this Chapter, the ID MERA for systems 
of anyons was also independently constructed by Konig and Bilgin (2010) on the torus. 
These authors provide proof of principle by computing ground state energies and two- 
point correlators for finite systems of Fibonacci anyons with x = [l? 1] (-^ = 2 in their 
notation), with errors in the energy on the order of a few percent. Once again we see 
that even for small values of the anyonic MERA is capable of providing an accurate 
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description of the low-energy behaviour of a system of interacting anyons. 

4.7 Some Notes on Implementation 
4.7.1 Block Structure 

Recall that in Sec. 3.5 we saw that an efficient way of storing Abelian symmetric tensors 
was as a number of blocks, and that operations on these tensors such as fusing and splitting 
legs amounted to nothing more than a rearrangement of these blocks. It is possible to 
implement a similar scheme for anyonic tensor networks, although with some important 
differences. 

As with the Abelian symmetric tensors of Chapter 3, we construct the central objects 
(e.g. M^^) of our tensors from a number of blocks. However, in contrast to the Abelian 
tensors, we no longer simply assemble these blocks into the final object. Instead, we give 
each block a unique identifying index, and introduce a map. This map then records the 
locations of the blocks, with each entry in the map corresponding to a unique labelling 
of all fusion trees (Fig. 4.22). Note that a given entry in the map may contain a list of 
more than one block, and each block listed is associated with a numeric multiplier. To 
reconstruct the tensor from the map and blocks, each entry in the map is assembled by 
summing over the relevant blocks, each multiplied by their associated numeric multiplier. 

This system may at first seem unwieldly. However, it has a number of advantages. First, 
F moves may be performed quickly and efficiently, with linear recombinations of the 
blocks of the tensor being performed simply by modifying the map. Second, braiding will 
permute the entries within individual blocks; efficiency gains may be made by deferring 
these permutations for as long as possible, letting them accumulate, and then determining 
and performing a single, cumulative operation on each block. By being able to perform 
F moves without accessing the contents of individual blocks, it becomes unnecessary 
to evaluate these deferred permutations when performing F moves (or indeed any other 
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Figure 4.22 - Map-based implementation of the central object Ma of an anyonic tensor: 
The diagram shows an operator M which acts on two adjacent sites of a ID lattice, where 
each site may carry charge I with degeneracy 1, or charge r with degeneracy 2. The associated 
tables show how the central object Ma of this particular operator is assembled in terms of its 
constituent blocks. If appropriate, a block may appear in more than one location in the map 
which describes Ma , and a location on the map may contain more than one block. 
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unitary operations depending only on the charges, such as the reversal of Frobenius-Schur 
indicators), and this can lead to greater computational efficiency. 

Further, the numerical factors associated with this braiding may also be introduced at 
the level of the map. This is typically an advantage when a particular charge labelling 
of the fusion trees is multiply degenerate, as is common at higher levels of the MERA. 
Each charge labelling is still only associated with a single block, and so applying these 
numerical factors at the level of the map requires less operations than applying them 
directly to the numerical content of the tensor. 

Finally, in Sec. 4.4.2 it was mentioned that it could be convenient to separate out the 
handling of normalisation factors associated with the diagrammatic isotopy convention. 
This takes place at two levels. First, when performing an operation such as "raising" 
an anyonic tensor to act on the fusion space of a larger number of sites, it is frequently 
possible to observe that factors coming from the introduction of vertices exactly cancel 
those arising due to the presence of loops; this is the reason why no factors of quantum 
dimensions da appear in (4.11). Second, when drawing an anyonic tensor with fusion trees, 
the vertex normalisation factors may be kept associated with the vertex. They therefore 
do not enter into M^^ except when M^'^ absorbs a vertex during the process of splitting 
or fusion of multi-indices. During splitting, when a vertex is absorbed, the associated 
factors of [dc/ {dadb)]^ are applied to the coefficients in the map, and not the entries of the 
tensor itself. During fusion, again the coefficients are applied only to the map, only this 
time they account for both the vertex itself and the loop which it makes with the central 
object of the tensor. 

With appropriate care, any effects of operations on the contents of individual blocks may 
be deferred until an operation is performed which by its nature must access the contents 
of the blocks, such as a singular value decomposition or a matrix multiplication. 

4.7.2 Precomputation 

If an algorithm employs repeated application of the same series of tensor manipulations, 
for example the repeated iterations of optimisation for the MERA, then many calculations 
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involved in these operations may be stored, and recycled on subsequent iterations. Exam- 
ples are the matrices generated by F moves and braids (Fig. 4.9), and the permutations 
of the numerical elements of each block which are generated by braiding. It is now that 
the ability to defer permutations of the elements of a tensor really comes into its own, as 
it is only necessary to store the cumulative operation for subsequent iterations, and not 
each individual step. 

It should be recognised that this Section describes only one particular scheme for the effi- 
cient implementation of anyonic tensors, and that this approach is by no means necessarily 
the only means of achieving this. However, for those interested in pursuing this approach 
further, additional discussion of the philosophy of precomputation [in this instance, as 
applied to U(l)-symmetric tensors] may be found in the Appendix of Singh et al. (2011), 
and Sec. 3.4 of this Thesis. 

4.7.3 Reminder: The Important Difference Between Fusion and 
Splitting Trees 

Given the implementation-oriented nature of this final Section of the Chapter, it seems 
appropriate to include a timely reminder of the difference between fusion and splitting 
trees. Recall that a tree assembled from fusion vertices describes the state of a system, 
IV'), and one assembled from splitting vertices describes a state in the dual space, 
where Hermitian conjugation is performed by vertical reflection of a tree diagram and 
complex conjugation of its coefficients. 

Note well that this same rule for Hermitian conjugation applies also to the vertex tensors 
A^2/3u ^iid N^"'^^, and thus fusing and splitting of legs takes place differently depending 
on whether it is acting on a fusing or a splitting tree. This important distinction will 
affect not only the arrangement of the charge blocks, but also of the entries within each 
block itself, and the correct implementation of these processes is one of the cornerstones 
for implementation of an anyonic tensor network. 
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Chapter 5 

Non-Abelian Symmetries of Spin Systems 



5.1 Unitary Braided Tensor Categories and Group 
Symmetries 

The formalism developed in Chapter 4 constitutes a methodology for performing tensor 
network simulations of any physical system which admits a description in terms of a 
UBTC. As mentioned in the introduction to that Chapter, it is also possible to associate 
a UBTC with a group Q, where the F moves are related to the 6-j symbols of the group, 
and the tensor Rf' is related to the choice of universal braid matrix (which describes the 
exchange properties of the irreps) . Note that the choice of universal braid matrix is not in 
general unique; in fact, we have already seen two systems with Z2 symmetry but different 
braiding. The first was the spin-0 formulation of the Ising model of Eq. (3.1), 

Rising = - 5^ (3.1) 

s 

where the Z2 symmetry is associated with a vr-radian rotation. If we write the charge 
labels of Z2 as (= +) and 1 (= — ), then the fusion rules may be written 

a X b ^ {a + b)\2 (5.1) 

(where I2 denotes that the addition is performed modulo 2), and the nonzero entries in 
are given by 

RU)U = 1- (5-2) 
163 
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The second example of a system with Z2 symmetry is any system of fermions, as discussed 
in Sec. 3.5.1, for example the fermionic formulation of the Ising model: 

^Ling = --^ E (c^^'^c(^+i) + ct(^+i)c(^) + ct('^)ct(^+i) + c(^+i)c(^) - 2(7ct(^)c(^) + g) . (5.3) 

s 

Here, the charges of correspond to parity, indicating the presence or absence of a 
fermion at a site s. The operators S^^^ and 6-^^ in the Hamiltonian are fermionic creation 
and annihilation operators, and the non-zero entries in i?"* for the associated UBTC are 
given by 

= (5.4) 
As Z2 is Abelian, the non-zero F moves in both examples are simply 



[There are no vertex indices u\, M2, M3, (see Fig. 4.1) in this expression, as all fusion 
products are non-degenerate.] 

For a non- Abelian group such as SU(2), both the F tensor and R^^ may be more compli- 
cated. The charges in SU(2) are the non- negative half- integers, with fusion rules 



a+6 

a X 6 — \ ^ c, (5.6) 

c=\a—h\ 



and we may write the F moves for SU(2) as 



a b e 



{Ft),f = (-l)^"^'+^+'V(2e + l)(2/ + l)<| } (5.7) 

c d f 
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a b e 

where < } denotes the 6-j symbol 

c d f 




A{a,b,e) A{e,c,d) A{b,c,f) A{a,f,d) 

{-lYiz + iy. 



[z — a — b — e)\{z — e — c — d)\{z — b — c — — a — f — d)\ 

1 



[a + b + c + d- z)\{a + e + c + f ~ z)\{b + e + d + f - z)\ 



A(a, 6, c) 



-a + b + c)\{a-b + c)!(a + b-c)\ 



(a + 6 + c+ 1)! 



(5.8) 
(5.9) 



with the sum running over all integer values of z such that the factorials are of non- 
negative numbers. When representing quantum mechanical spin, the half-integer charges 
are fermionic, and thus the tensor is given by 



for any combination of a, b, and c permitted by the fusion rules (5.6). 



(5.10) 



Given the F tensor and the particle exchange tensor Rf', we may apply the UBTC formal- 
ism of Chapter 4 to any quantum mechanical system which exhibits a group symmetry. 
For Abelian symmetry groups, much of this machinery is redundant and we may prefer 
the simpler approach outlined in Chapter 3. For non-Abelian symmetries, however, this 
provides a useful means of exploiting those symmetries. 

Expressions (5.1), (5.2), and (5.4)-(5.10) in this Section are adapted from expressions 
found in Chapter 5 of Bonderson (2007). 



5.2 Fermions revisited 



Before presenting the application of the formalism of Chapter 4 to an example of a system 
exhibiting a non-Abelian symmetry, we will first consider its application to a system of 
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fermions. In Sec. 3.5.1 we saw that adding fermionic statistics to a symmetric tensor 
network algorithm involved introducing extra factors of —1 into permutation operations, 
and also into the multiplications used to perform tensor contraction. This process is 
counter-intuitive and an approach would be preferable in which factors associated with 
particle exchange arise only during index permutation. This may achieved by representing 
fermionic systems using the F and tensors of Eqs. (5.5) and (5.4) in the UBTC 
formalism of Chapter 4. However, in truth, the full machinery of the UBTC tensor network 
formalism is not necessary. Instead, it suffices to simply assign a vertical orientation to 
each leg and to contract pairs of upgoing legs using the expanded multiplicity tensor N]^^, 
and pairs of downgoing legs using its Hermitian conjugate N^'^^ , as per Sec. 4.4.2. When 
contracting two tensors together, the counterclockwise ordering of indices on one of these 
tensors is now the opposite of that on the other, e.g. 

u 



u I 

9 

w I 




r" = b[ ]c , (5.11) 

w ' 

/ 

and there are no longer any concealed particle exchanges within the equivalent of Fig. 3.8(3), 

(5.12) 



Of" 

w I 




which is now just a simple matrix multiplication. 



5.3 The ID Spin-| Heisenberg (Anti)Ferromagnet 

The procedure for studying systems with non-Abelian symmetries using a UBTC-based 
tensor network is largely the same as that described in Chapter 4, with only a couple of 
minor differences. First, until now I have implicitly employed a convention where every 
fusion tree has a total charge of I, and represents a single state. It is also possible to use 
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a fusion tree with a total charge a 7^ I to represent a single state, where a total charge 
of a is taken to imply the existence of an ancillary system not explicitly considered, and 
having charge a (for example, for a system of anyons on a finite disc this ancillary charge 
may live on the boundary). However, in this Section it is instead preferable to use a 
fusion tree with a total charge of a to represent a subspace of the Hilbert space having 
dimension da} Second, if we are considering systems with only fermionic and/or bosonic 
statistics, it is not necessary to specify the orientation of a braid. Thus we may denote 
particle exchange simply by line crossings, 




(5.13) 



«3 «3 «3 



As examples of systems with SU(2) invariance, we shall consider the ID spin-| Heisenberg 
antiferromagnet and ferromagnet, with periodic boundary conditions. The former model 
exhibits a nearest-neighbour Hamiltonian which favours neighbouring pairs of particles 
occupying the singlet state, with total spin 0, 



"■AFM ~ 





A 

1 1 

2 2 



(5.14) 



^There is a subtlety here, in that for a system with SU(2) symmetry, we might choose to define the 
Hilbert space such that an orthonormal basis is given by eigenstates of a complete set of SU(2)-symmetric 
commuting operators. Under this choice, a fusion tree always corresponds to a single state, regardless 
of total charge. More commonly, however, for systems exhibiting a non-Abelian symmetry we define 
the Hilbert space with respect to a complete set of commuting operators on the microscopic degrees of 
freedom of the system. These measurements are not necessarily SU(2)-symmetric, e.g. measurement of 
spin in the z basis for a Heisenberg spin chain. In the resulting basis of this example, a pair of spin-1/2 
fermions can have total spin in precisely one way, or total spin 1 in three orthogonal ways, corresponding 
to the SU(2)-symmetry-breaking measurement of z-axis spin permitting resolution of a three-dimensional 
subspace for the spin-1 triplet [which has a total SU(2) "charge" of 1, with di =3]. 
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whereas the latter favours occupation of the triplet sector, with spin 1, 



h 



FM 



On two sites, the total Hamiltonians may be written as 



(5.15) 



Hafm 





A 





(5.16) 



1 

A 





(5.17) 



and the diagrams 



\y and \y 

'o 'l 



(5.18) 



correspond to eigensubspaccs of these Hamiltonians, with the former representing the 
spin-0 singlet state [total SU(2) "charge" 0; cIq = 1, hence a state] and the latter rep- 
resenting the 3-dimensional spin-1 triplet sector [total SU(2) "charge" 1; di = 3, hence 
a 3-dimensional subspace]. For the given antiferromagnetic Hamiltonian these diagrams 
have eigenvalues —2 and respectively, whereas for the ferromagnetic Hamiltonian these 
eigenvalues are reversed. 
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This concludes a very simple demonstration of the application of UBTC tensor networks 
for the exploitation of non-Abelian symmetries. A fuller treatment will be provided in 
Singh and Vidal (in preparation). 
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Chapter 6 



Summary and Outlook 



In this Thesis, we have seen how symmetries — both spatial and internal — may be exploited 
in tensor network algorithms. We began in Chapter 2 with the exploitation of scale 
invariance, constructing an Ansatz which naturally reflects the entanglement structure 
present in quantum critical systems. Applying this Ansatz to the critical Ising and three- 
state Potts models in ID, we were able to extract most of the conformal data of the CFTs 
which describe the continuum limit of these theories, namely the scaling dimensions, 
central charge, and the coefficients of the operator product expansion. 

In Chapter 3 we described the mathematical background behind internal symmetries of 
lattice models, and developed basic techniques to exploit Abelian symmetries in tensor 
network algorithms. These techniques enabled us to address specific symmetry sectors 
of models such as the XX and Heisenberg models, and to simulate these systems at a 
substantially reduced computational cost. In their own rights, these techniques therefore 
substantially increase the power of tensor network algorithms as tools for the study of 
symmetric systems. 

In Chapter 4, we turned our attention to systems of anyons. The simulation of anyonic 
systems is acknowledged to be challenging: Like fermions, the study of anyons suffers 
from the sign problem,^ and the problem of developing a general algorithm to compute 

^The sign problem may be paraphrased as "For the class of Hamiltonians obtained by taking a bosonic 
system whose ground state energy may be computed in polynomial time and introducing fermionic particle 
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the ground state of arbitrary fermionic systems in polynomial time is therefore known 
to be at least NP-hard (Troycr and Wiese, 2005). Nevertheless, by means of a non- 
trivial generalisation of the techniques introduced in Chapter 3, we were able to develop 
a formalism of tensor networks for anyons, allowing us to compute a close approximation 
to the ground state of an anyonic system in polynomial time, provided the entanglement 
structure of that ground state may be effectively represented by an appropriate tensor 
network algorithm (see the discussion on entanglement in Sec. 2.1.2). 

Finally, in Chapter 5 it was seen that the formalism developed for anyons in Chapter 4 
may also be applied to exploit non-Abelian symmetries of spin systems, such as the SU(2) 
symmetry of the Heisenberg (anti)ferromagnet, or indeed to exploit Abelian symmetries in 
the presence of possibly non-trivial exchange statistics. In fact, the formalism of Chapter 4 
may be used to study systems of bosons, fermions, Abelian and non-Abelian anyons, and 
to exploit the presence of Abelian and non-Abelian internal symmetries of the Hamiltonian 
in any of these systems.'^ 

The exploitation of symmetries is a powerful tool, vastly increasing the reach and power 
of tensor network algorithms as a condensed matter technique. Of all the developments 
described above, perhaps the most exciting is the extension of tensor network algorithms 
to anyons, opening the door to the study of a vast array of condensed matter systems, 
many of which have never been studied before. Many of the questions to be asked are of 
great significance — for example, consider the Fibonacci anyons studied in Sec. 4.5.4. They 
can implement universal quantum computation through braiding alone, are believed to 



exchange statistics, does there exist an algorithm similarly capable of computing the ground state energy 
of these fermionic Hamiltonians in polynomial time?" This question, which is stated more formally 
Troycr and Wiese (2005), remains unanswered, and as demonstrated by Troycr and Wiese, is in fact 
NP-hard. Specific solutions are known for many problems and problem groups, but there exist many 
other such fermionic systems whose ground state may at present only be computed for a cost exponential 
in the system size, even though their bosonised equivalents may be solved in polynomial time. There 
arc of course bosonic systems (such as frustrated systems, and spin glasses) which are also exponentially 
hard, and their fermionic counterparts tend to be exponentially hard as well. It is the presence of the 
middle ground, where fermions are "harder" than bosons, to which the term "the sign problem" is usually 
applied. 

^With minimal modification, the formalism may even be applied to systems where particle exchange 
is not possible, and the system is a ID open chain described by a unitary tensor category admitting 
no solutions to the hexagon equation; however, the author is as yet unaware of any interesting physical 
models of this form. 
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appear as quasiparticles in the v = 12/5 fractional quantum Hall state, and yet we are only 
beginning to understand their phase diagrams under even the simplest of interactions (e.g. 
Trebst et al., 2008a). Anyonic tensor networks are a powerful tool for asking fundamental 
questions about such systems, and could be of vital importance to the coming quantum 
revolution in information processing. They also provide an unrivalled opportunity to gain 
insight into this fascinating and comparatively little-understood area of condensed matter 
physics. 



This is an exciting time to be working on anyons! 
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